## [1] "2019-01-27-17-08-53"
## [1] "/home/peterslab/Alex Alleman/Statewide Microbiome Analysis/Statewide analysis"
set.seed(8765)
library(ggplot2)
library(data.table)
library(vegan)
library(ggvegan)
library(dplyr)
library(scales)
library(grid)
library(reshape2)
library(phyloseq)
library(ggpubr)
library(RColorBrewer)
library(ape)
library(grid)
library(knitr)
library(igraph)
library(Matrix)
library(ggnetwork)
library(intergraph)
library(parallel)
library(tinytex)
farm_col<-(c("#8c510a", "#d8b365", "#f6e8c3", "#f5f5f5", "#c7eae5", "#5ab4ac", "#01665e"))
farm_col_dark<-brewer.pal(7, "Dark2")
farm_col_paired<-brewer.pal(10, "Paired")
OTU_nifH<- read.delim(
"~/Alex Alleman/Statewide Microbiome Analysis/Statewide analysis/nifH_OTUallspring.csv",
row.names = 1)
head(OTU_nifH)[,1:10]
## JZ017 JZ018 JZ019 JZ020 JZ021 JZ022 JZ023 JZ024 JZ025 JZ026
## OTU1 21 32923 48199 22 117 103 147 4054 5597 246
## OTU2 35281 5 15 1075 101 1335 26249 376 1705 256
## OTU3 11 14 13 49 20 48 9 8 16 9
## OTU4 80 4 14 2576 144 136 150 129 29299 77807
## OTU5 6 2 7 146 749 526 432 1115 13256 0
## OTU6 19317 11274 799 748 1532 348 32 21 28 2517
tax_nifH<- read.delim(
"~/Alex Alleman/Statewide Microbiome Analysis/Statewide analysis/nifH_OTU_ids_2016_fixed1.txt",
row.names = 1)
head(tax_nifH)[,1:8]
## kingdom phylum class order
## OTU1 k__bacteria p__firmicutes c__bacilli o__bacillales
## OTU2 k__bacteria p__firmicutes c__bacilli o__bacillales
## OTU3 k__bacteria p__proteobacteria c__alphaproteobacteria o__rhizobiales
## OTU4 k__bacteria p__firmicutes c__bacilli o__bacillales
## OTU5 k__bacteria p__proteobacteria c__alphaproteobacteria o__rhizobiales
## OTU6 k__bacteria p__proteobacteria c__alphaproteobacteria o__rhizobiales
## family genus
## OTU1 f__paenibacillaceae g__paenibacillus
## OTU2 f__paenibacillaceae g__paenibacillus
## OTU3 f__rhizobiaceae g__rhizobium
## OTU4 f__paenibacillaceae g__paenibacillus
## OTU5 f__bradyrhizobiaceae g__bradyrhizobium
## OTU6 f__bradyrhizobiaceae g__rhodopseudomonas
## species strain
## OTU1 s__paenibacillus graminis paenibacillus graminis
## OTU2 s__paenibacillus polymyxa paenibacillus polymyxa
## OTU3 s__rhizobium leguminosarum rhizobium leguminosarum
## OTU4 s__paenibacillus terrae paenibacillus terrae hpl_003
## OTU5 s__bradyrhizobium sp. bradyrhizobium sp. stm3074
## OTU6 s__rhodopseudomonas palustris rhodopseudomonas palustris
meta<- read.delim(
"~/Alex Alleman/Statewide Microbiome Analysis/Statewide analysis/all_metadata_summer.csv", row.names = 1, na.strings = "NA", sep = ",",
colClasses = c(rep('factor', 7), rep('numeric', 3), rep('factor', 4), 'numeric',
rep('factor', 3), rep('numeric', 28) ) )
head(meta)[,1:5]
## Site ARC Season Sample_dates Pea_variety
## JZ032 Kalispell NWARC Summer 2016-summer Delta
## JZ031 Kalispell NWARC Summer 2016-summer CDC Saffron
## JZ030 Kalispell NWARC Summer 2016-summer AC Earlystar
## JZ034 Kalispell NWARC Summer 2016-summer Majoret
## JZ033 Kalispell NWARC Summer 2016-summer DS Admiral
## JZ035 Kalispell NWARC Summer 2016-summer Navarro
meta2 <- meta[-c(19:48),]
sapply(meta2, class)
## Site ARC Season Sample_dates
## "factor" "factor" "factor" "factor"
## Pea_variety Plot season_precip irrgation
## "factor" "factor" "numeric" "numeric"
## total_precip_irr sample_depth Date Tillage
## "numeric" "factor" "factor" "factor"
## prev_crop grain_yield elevation lat
## "factor" "numeric" "factor" "factor"
## lon Organic_Matter Moisture_Content Nitrate_Nitrite
## "factor" "numeric" "numeric" "numeric"
## Ammonia Av_Phosphorus Av_Potassium Sulfate_Sulfur
## "numeric" "numeric" "numeric" "numeric"
## pH Boron Arsenic Barium
## "numeric" "numeric" "numeric" "numeric"
## Cadmium Calcium Chromium Cobalt
## "numeric" "numeric" "numeric" "numeric"
## Copper Iron Lead Magnesium
## "numeric" "numeric" "numeric" "numeric"
## Manganese Molybdenum Nickel Phosphorus
## "numeric" "numeric" "numeric" "numeric"
## Potassium Sodium Sulfur Vanadium
## "numeric" "numeric" "numeric" "numeric"
## Zinc
## "numeric"
OTU_nifH_m<-as.matrix(OTU_nifH)
tax_nifH_m<-as.matrix(tax_nifH)
meta_m<-as.matrix(meta2)
class(OTU_nifH_m)
## [1] "matrix"
class(tax_nifH_m)
## [1] "matrix"
class(meta_m)
## [1] "matrix"
OTUnifH = otu_table(OTU_nifH_m, taxa_are_rows = TRUE)
TAXnifH = tax_table(tax_nifH_m)
physeq_nifH = phyloseq(OTUnifH, TAXnifH)
physeq_nifH
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 8821 taxa and 101 samples ]
## tax_table() Taxonomy Table: [ 8821 taxa by 8 taxonomic ranks ]
meta_phy <- sample_data(meta2)
sample_names(meta_phy)
## [1] "JZ032" "JZ031" "JZ030" "JZ034" "JZ033" "JZ035" "JZ040" "JZ046"
## [9] "JZ042" "JZ044" "JZ038" "JZ036" "JZ041" "JZ047" "JZ043" "JZ045"
## [17] "JZ039" "JZ037" "JZ081" "JZ078" "JZ082" "JZ079" "JZ080" "JZ083"
## [25] "JZ105" "JZ107" "JZ103" "JZ102" "JZ106" "JZ104" "JZ084" "JZ085"
## [33] "JZ086" "JZ087" "JZ088" "JZ089" "JZ090" "JZ091" "JZ092" "JZ093"
## [41] "JZ094" "JZ095" "JZ096" "JZ097" "JZ098" "JZ099" "JZ100" "JZ101"
## [49] "JZ112" "JZ109" "JZ110" "JZ111" "JZ108" "JZ113"
physeq_nifH<-merge_phyloseq(physeq_nifH, meta_phy)
physeq_nifH
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 8821 taxa and 54 samples ]
## sample_data() Sample Data: [ 54 samples by 45 sample variables ]
## tax_table() Taxonomy Table: [ 8821 taxa by 8 taxonomic ranks ]
Source of trim protocol http://evomics.org/wp-content/uploads/2016/01/phyloseq-Lab-01-Answers.html#taxa-total-counts-histogram
tdt_nifH = data.table(tax_table(physeq_nifH),
TotalCounts = taxa_sums(physeq_nifH),
OTU = taxa_names(physeq_nifH))
ggplot(tdt_nifH, aes(TotalCounts)) +
geom_histogram() +
ggtitle("Histogram of Total Counts nifH")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# How many OTUS have low count (Rare)?
tdt_nifH[(TotalCounts <= 0), .N]#zero count
## [1] 1691
tdt_nifH[(TotalCounts <= 1), .N]#single count
## [1] 2337
tdt_nifH[(TotalCounts <= 2), .N]#double count
## [1] 2778
# taxa cumulative sum
taxcumsum_nifH = tdt_nifH[, .N, by = TotalCounts]
setkey(taxcumsum_nifH, TotalCounts)
taxcumsum_nifH[, CumSum := cumsum(N)]
# Define the plot
pCumSum_nifH = ggplot(taxcumsum_nifH, aes(TotalCounts, CumSum)) +
geom_point() +
xlab("Filtering Threshold, Minimum Total Counts") +
ylab("OTUs Filtered") +
ggtitle("OTUs that would be filtered vs. the minimum count threshold (nifH)")
pCumSum_nifH
pCumSum_nifH + xlim(0, 100)
## Warning: Removed 856 rows containing missing values (geom_point).
mdt_nifH = fast_melt(physeq_nifH)
prevdtnifH = mdt_nifH[, list(Prevalence = sum(count > 0),
TotalCounts = sum(count)),
by = TaxaID]
ggplot(prevdtnifH, aes(Prevalence)) +
geom_histogram() +
ggtitle("Histogram of Taxa Prevalence nifH")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# How many OTUS have low prevelance (Rare)?
prevdtnifH[(Prevalence <= 0), .N]#zero
## [1] 1691
prevdtnifH[(Prevalence <= 1), .N]#single
## [1] 3277
prevdtnifH[(Prevalence <= 2), .N]#double
## [1] 4476
prevcumsumnifH = prevdtnifH[, .N, by = Prevalence]
setkey(prevcumsumnifH, Prevalence)
prevcumsumnifH[, CumSum := cumsum(N)]
pPrevCumSumnifH = ggplot(prevcumsumnifH, aes(Prevalence, CumSum)) +
geom_point() +
xlab("Filtering Threshold, Prevalence") +
ylab("OTUs Filtered") +
ggtitle("OTUs that would be filtered vs. the minimum count threshold")
pPrevCumSumnifH
ggplot(prevdtnifH, aes(Prevalence, TotalCounts)) +
geom_point(size = 2, alpha = 0.5) +
scale_y_log10()
*** ###Trimming ####Remove less than doublets in data and prevlant in 5% of the sample ####First transform to realtive abundance
physeq_nifH_trim = filter_taxa(physeq_nifH, function(x) sum(x > 2) > (0.05*length(x)),
TRUE)
physeq_nifH_trim
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 2107 taxa and 54 samples ]
## sample_data() Sample Data: [ 54 samples by 45 sample variables ]
## tax_table() Taxonomy Table: [ 2107 taxa by 8 taxonomic ranks ]
wh0_n = genefilter_sample(physeq_nifH_trim, filterfun_sample(function(x) x > 5),
A=0.1*nsamples(physeq_nifH_trim))
physeq_nifH_ord = prune_taxa(wh0_n, physeq_nifH_trim)
physeq_nifH_ord
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 804 taxa and 54 samples ]
## sample_data() Sample Data: [ 54 samples by 45 sample variables ]
## tax_table() Taxonomy Table: [ 804 taxa by 8 taxonomic ranks ]
physeq_nifH_ord = transform_sample_counts(physeq_nifH_ord, function(x) 1E6 * x/sum(x))
physeq_nifH_ord
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 804 taxa and 54 samples ]
## sample_data() Sample Data: [ 54 samples by 45 sample variables ]
## tax_table() Taxonomy Table: [ 804 taxa by 8 taxonomic ranks ]
physeq_nifH_ord_2_phylum <- tax_glom(physeq_nifH_ord, "phylum")
plot_bar(physeq_nifH_ord_2_phylum, x = "Site", fill = "phylum")+
geom_bar(aes(fill = phylum), color = "black", stat = "identity", position = "stack",
show.legend = TRUE)+
scale_fill_manual(name = "Phylum",
values=c("#7fc97f", "#beaed4", "#fdc086", "#ffff99", "#386cb0",
"#f0027f",'#e31a1c'),
#breaks=c("p_actinobacteria", "p_bacteroidetes", "p_cyanobacteria",
#"p_firmicutes", "p_proteobacteria", "p_verrucomicrobia"),
labels = c("Actinobacteria", "Bacteroidetes", "Cyanobacteria",
"Firmicutes", "Proteobacteria", "Verrucomicrobia"),
guide = guide_legend(reverse = TRUE)
)+
ggtitle("nifH Phylum Relative Abundance by Site")+
ylab("Relative Abundance (%)")+
theme(axis.text.x = element_text(angle = 90, hjust = 1), axis.text.y = element_blank(),
panel.background = element_blank())
physeq_nifH_ord_2_phylum <- tax_glom(physeq_nifH_ord, "phylum")
plot_bar(physeq_nifH_ord_2_phylum, x = "Site", fill = "phylum")+
geom_bar(aes(fill = phylum), stat = "identity", position = "stack", show.legend = TRUE)+
scale_fill_manual(name = "Phylum",
values=c("#7fc97f", "#beaed4", "#fdc086", "#ffff99",
"#386cb0", "#f0027f"),
#breaks=c("p_actinobacteria", "p_bacteroidetes", "p_cyanobacteria",
#"p_firmicutes", "p_proteobacteria", "p_verrucomicrobia"),
labels = c("Actinobacteria", "Bacteroidetes", "Cyanobacteria",
"Firmicutes", "Proteobacteria", "Verrucomicrobia"),
guide = guide_legend(reverse = TRUE)
)+
ggtitle("nifH Phylum Relative Abundance by Site")+
ylab("Relative Abundance (%)")+
theme(axis.text.x = element_text(angle = 90, hjust = 1), axis.text.y = element_blank(),
panel.background = element_blank())
physeq_nifH_ord_1 = transform_sample_counts(physeq_nifH_ord, function(x) x / sum(x) )
physeq_nifH_ord_phylum <- tax_glom(physeq_nifH_ord_1, "phylum")
data_nifH_phylum <- psmelt(physeq_nifH_ord_phylum)
data_nifH_phylum$phylum<-as.character(data_nifH_phylum$phylum)
data_nifH_phylum$phylum[data_nifH_phylum$Abundance<0.05]<-"<5% abdund"
#list new phylums
unique(data_nifH_phylum$phylum)
## [1] "p__proteobacteria" "p__firmicutes" "p__actinobacteria"
## [4] "<5% abdund"
ggplot(data = data_nifH_phylum, aes(x = Site, y = Abundance, fill = phylum))+
geom_bar(aes(fill = phylum), stat = "identity", position = "stack", show.legend = TRUE)+
scale_fill_manual(name = "Phylum",
values=c( "#33a02c", "#b2df8a","#a6cee3", "#1f78b4"),
breaks=c("p__proteobacteria", "p__firmicutes",
"p__actinobacteria","<5% abdund" ),
labels = c( "Proteobacteria", "Firmicutes",
"Actinobacteria", "<5% Abundance"),
guide = guide_legend(reverse = TRUE)
)+
ggtitle("nifH Phylum Relative Abundance by Site")+
ylab("Relative Abundance (%)")+
scale_x_discrete(labels = c("Conrad", "Corvallis", "Huntley Dryland",
"Huntley Irrigated", "Kalispell", "Mocassin",
"Sidney Dryland", "Sidney Irrigated", "Richland"))+
theme(axis.text.x = element_text(angle = 90, hjust = 1),
axis.text.y = element_blank(), panel.background = element_blank())
## png
## 2
ggplot(data = data_nifH_phylum, aes(x = Site, y = Abundance, fill = phylum))+
geom_bar(aes(fill = phylum), stat = "identity", position = "stack", show.legend = TRUE)+
scale_fill_manual(name = "Phylum",
values=c( "#33a02c", "#b2df8a","#a6cee3", "#ff7f00"),
breaks=c("p__proteobacteria", "p__firmicutes",
"p__actinobacteria","<5% abdund" ),
labels = c( "Proteobacteria", "Firmicutes",
"Actinobacteria", "<5% Abundance"),
guide = guide_legend(reverse = TRUE)
)+
ggtitle("nifH Phylum Relative Abundance by Site")+
ylab("Relative Abundance")+
scale_x_discrete(labels = c("Conrad", "Corvallis",
"Huntley Dryland", "Huntley Irrigated", "Kalispell",
"Mocassin", "Sidney Dryland", "Sidney Irrigated", "Richland"))+
theme(axis.text.x = element_text(angle = 90, hjust = 1),
axis.text.y = element_blank(), panel.background = element_blank())
physeq_nifH_ord_1 = transform_sample_counts(physeq_nifH_ord, function(x) x / sum(x) )
physeq_nifH_ord_genus <- tax_glom(physeq_nifH_ord_1, "genus")
data_nifH_genus <- psmelt(physeq_nifH_ord_genus)
data_nifH_genus$genus<-as.character(data_nifH_genus$genus)
data_nifH_genus$genus[data_nifH_genus$Abundance<0.1]<-"<10% abdund"
unique(data_nifH_genus$genus)
## [1] "g__paenibacillus" "g__rhizobium" "g__geobacter"
## [4] "g__rhodopseudomonas" "g__bradyrhizobium" "g__gluconacetobacter"
## [7] "g__dechloromonas" "g__agrobacterium" "g__ideonella"
## [10] "g__azospirillum" "g__arthrobacter" "<10% abdund"
ggplot(data = data_nifH_genus, aes(x = Site, y = Abundance, fill = genus))+
geom_bar(aes(fill = genus), stat = "identity", position = "stack", show.legend = TRUE)+
scale_fill_manual(name = "Genus",
values=c('#fdbf6f','#a6cee3','#b2df8a','#fb9a99','#e31a1c','#ff7f00', '#cab2d6','#6a3d9a','#ffff99','#80b1d3','#1f78b4', '#33a02c'),
breaks=c("g__paenibacillus", "g__rhizobium", "g__geobacter",
"g__rhodopseudomonas","g__bradyrhizobium",
"g__gluconacetobacter","g__dechloromonas",
"g__agrobacterium", "g__ideonella",
"g__azospirillum","g__arthrobacter", "<10% abdund"),
labels=c("Paenibacillus", "Rhizobium", "Geobacter", "Rhodopseudomonas", "Bradyrhizobium", "Gluconacetobacter", "Dechloromonas", "Agrobacterium", "Ideonella", "Azospirillum","Arthrobacter",
"<10% abdund"),
guide = guide_legend(reverse = TRUE)
)+
ggtitle("nifH Genus Relative Abundance by Site")+
ylab("Relative Abundance")+
scale_x_discrete(labels = c("Conrad", "Corvallis", "Huntley Dryland",
"Huntley Irrigated", "Kalispell", "Mocassin",
"Sidney Dryland", "Sidney Irrigated", "Richland"))+
theme(axis.text.x = element_text(angle = 90, hjust = 1), axis.text.y = element_blank(),
panel.background = element_blank())
## png
## 2
plot_richness(physeq_nifH_trim)
plot_richness(physeq_nifH_trim, measures = c("Observed","Chao1", "Shannon"), color = "Site")
rich_nifH<-estimate_richness(physeq_nifH_trim, split = TRUE)
#First merge data sets with meta2
meta2$sample_names<-rownames(meta2)
rich_nifH$sample_names<-rownames(rich_nifH)
meta_nifH<-merge(meta2, rich_nifH, by = "sample_names")
rownames(meta_nifH)<-meta_nifH$sample_names
meta_nifH<-meta_nifH[,-1]
head(meta_nifH)[,49:54]
## ACE se.ACE Shannon Simpson InvSimpson Fisher
## JZ030 1061.8638 16.02588 3.840892 0.9370300 15.880569 130.0829
## JZ031 1015.7345 15.87722 3.786340 0.9379526 16.116700 120.5734
## JZ032 985.3458 15.23833 3.241100 0.8739992 7.936459 116.9659
## JZ033 1077.3023 16.24837 3.968293 0.9463960 18.655339 133.1658
## JZ034 1044.4015 15.97121 3.962209 0.9532236 21.378282 129.1375
## JZ035 1123.6740 16.58229 3.653639 0.9231705 13.015828 133.4697
#use ggpubr for plot
nifH_Observ<-ggboxplot(meta_nifH, x = "Site", y = "Observed",
rug = TRUE,
fill = "Site", xlab = " ",
palette = farm_col_paired)+
rremove("x.text")
nifH_Shannon<-ggboxplot(meta_nifH, x = "Site", y = "Shannon",
rug = TRUE,
fill = "Site", xlab = " ",
palette = farm_col_paired)+
rremove("x.text")
nifH_Chao<- ggboxplot(meta_nifH, x = "Site", y = "Chao1",
rug = TRUE,
fill = "Site", xlab = " ",
palette = farm_col_paired)+
rremove("x.text")
nifH_InvSim<- ggboxplot(meta_nifH, x = "Site", y = "InvSimpson",
rug = TRUE,
fill = "Site", xlab = " ",
palette = farm_col_paired)+
rremove("x.text")
alpha_nifH_fig<-ggarrange(nifH_Observ, nifH_Shannon, nifH_Chao, nifH_InvSim, ncol = 2, nrow = 2, common.legend = TRUE, legend = "bottom")
annotate_figure(alpha_nifH_fig, top = text_grob("Alpha Diversity of nifH by Site", size = 20))
Publish to .tiff
## png
## 2
Source: http://biology.kenyon.edu/courses/biol229/diversity.pdf
#use ggpubr for plot
ggboxplot(meta_nifH, x = "Site", y = "Shannon",
add = "mean", rug = TRUE,
fill = "Site",
title = "nifH Shannon Diversity by Site", palette = farm_col_paired, legend = "right")+
stat_compare_means(label.y = 4.3, p.adjust.method = "bonferroni")+
rotate_x_text()
## Warning: Ignoring unknown parameters: p.adjust.method
https://rpubs.com/dillmcfarlan/R_microbiotaSOP
#Create 2x2 plot environment so that we can see all 4 metrics at once.
par(mfrow = c(2, 2))
#Then plot each metric.
hist(rich_nifH$Shannon, main="Shannon diversity", xlab="", breaks=10)
hist(rich_nifH$InvSimpson, main="Inverse Simpson diversity", xlab="", breaks=15)
hist(rich_nifH$Chao1, main="Chao richness", xlab="", breaks=15)
hist(rich_nifH$ACE, main="ACE richness", xlab="", breaks=15)
shapiro.test(rich_nifH$Shannon)
##
## Shapiro-Wilk normality test
##
## data: rich_nifH$Shannon
## W = 0.97136, p-value = 0.2211
shapiro.test(rich_nifH$InvSimpson)
##
## Shapiro-Wilk normality test
##
## data: rich_nifH$InvSimpson
## W = 0.8798, p-value = 6.121e-05
shapiro.test(rich_nifH$Chao1)
##
## Shapiro-Wilk normality test
##
## data: rich_nifH$Chao1
## W = 0.94469, p-value = 0.01469
shapiro.test(rich_nifH$ACE)
##
## Shapiro-Wilk normality test
##
## data: rich_nifH$ACE
## W = 0.93376, p-value = 0.005166
aov_shannon_site_nifH <- aov(Shannon ~ Site, meta_nifH)
summary(aov_shannon_site_nifH)
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 8 15.624 1.9531 15.43 1.28e-10 ***
## Residuals 45 5.696 0.1266
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(aov_shannon_site_nifH)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Shannon ~ Site, data = meta_nifH)
##
## $Site
## diff lwr upr
## Corvallis-Conrad 0.29300868 -0.376019931 0.962037285
## Huntley Dryland -Conrad -0.50522512 -1.174253731 0.163803484
## Huntley Irrigated-Conrad 0.86418474 0.195156134 1.533213349
## Kalispell-Conrad 0.66189448 -0.007134125 1.330923091
## Mocassin-Conrad -0.44378350 -1.112812113 0.225245103
## Richland-Conrad -0.68161691 -1.350645518 -0.012588302
## Sidney Dryland-Conrad -0.66014194 -1.329170551 0.008886665
## Sidney Irrigated-Conrad 0.03893720 -0.630091409 0.707965807
## Huntley Dryland -Corvallis -0.79823380 -1.467262408 -0.129205192
## Huntley Irrigated-Corvallis 0.57117606 -0.097852543 1.240204673
## Kalispell-Corvallis 0.36888581 -0.300142802 1.037914414
## Mocassin-Corvallis -0.73679218 -1.405820790 -0.067763574
## Richland-Corvallis -0.97462559 -1.643654195 -0.305596979
## Sidney Dryland-Corvallis -0.95315062 -1.622179227 -0.284122012
## Sidney Irrigated-Corvallis -0.25407148 -0.923100086 0.414957130
## Huntley Irrigated-Huntley Dryland 1.36940986 0.700381257 2.038438473
## Kalispell-Huntley Dryland 1.16711961 0.498090998 1.836148214
## Mocassin-Huntley Dryland 0.06144162 -0.607586989 0.730470227
## Richland-Huntley Dryland -0.17639179 -0.845420395 0.492636821
## Sidney Dryland-Huntley Dryland -0.15491682 -0.823945427 0.514111789
## Sidney Irrigated-Huntley Dryland 0.54416232 -0.124866286 1.213190930
## Kalispell-Huntley Irrigated -0.20229026 -0.871318867 0.466738349
## Mocassin-Huntley Irrigated -1.30796825 -1.976996854 -0.638939638
## Richland-Huntley Irrigated -1.54580165 -2.214830260 -0.876773044
## Sidney Dryland-Huntley Irrigated -1.52432668 -2.193355292 -0.855298076
## Sidney Irrigated-Huntley Irrigated -0.82524754 -1.494276151 -0.156218935
## Mocassin-Kalispell -1.10567799 -1.774706596 -0.436649380
## Richland-Kalispell -1.34351139 -2.012540001 -0.674482785
## Sidney Dryland-Kalispell -1.32203643 -1.991065033 -0.653007818
## Sidney Irrigated-Kalispell -0.62295728 -1.291985892 0.046071324
## Richland-Mocassin -0.23783341 -0.906862013 0.431195203
## Sidney Dryland-Mocassin -0.21635844 -0.885387046 0.452670170
## Sidney Irrigated-Mocassin 0.48272070 -0.186307904 1.151749312
## Sidney Dryland-Richland 0.02147497 -0.647553640 0.690503576
## Sidney Irrigated-Richland 0.72055411 0.051525501 1.389582717
## Sidney Irrigated-Sidney Dryland 0.69907914 0.030050533 1.368107749
## p adj
## Corvallis-Conrad 0.8815206
## Huntley Dryland -Conrad 0.2786584
## Huntley Irrigated-Conrad 0.0035559
## Kalispell-Conrad 0.0545026
## Mocassin-Conrad 0.4476016
## Richland-Conrad 0.0428516
## Sidney Dryland-Conrad 0.0556618
## Sidney Irrigated-Conrad 0.9999999
## Huntley Dryland -Corvallis 0.0091627
## Huntley Irrigated-Corvallis 0.1498601
## Kalispell-Corvallis 0.6844416
## Mocassin-Corvallis 0.0211504
## Richland-Corvallis 0.0006683
## Sidney Dryland-Corvallis 0.0009315
## Sidney Irrigated-Corvallis 0.9435800
## Huntley Irrigated-Huntley Dryland 0.0000011
## Kalispell-Huntley Dryland 0.0000308
## Mocassin-Huntley Dryland 0.9999977
## Richland-Huntley Dryland 0.9940451
## Sidney Dryland-Huntley Dryland 0.9975503
## Sidney Irrigated-Huntley Dryland 0.1957237
## Kalispell-Huntley Irrigated 0.9854512
## Mocassin-Huntley Irrigated 0.0000030
## Richland-Huntley Irrigated 0.0000001
## Sidney Dryland-Huntley Irrigated 0.0000001
## Sidney Irrigated-Huntley Irrigated 0.0062523
## Mocassin-Kalispell 0.0000835
## Richland-Kalispell 0.0000017
## Sidney Dryland-Kalispell 0.0000024
## Sidney Irrigated-Kalispell 0.0858435
## Richland-Mocassin 0.9611572
## Sidney Dryland-Mocassin 0.9778887
## Sidney Irrigated-Mocassin 0.3356109
## Sidney Dryland-Richland 1.0000000
## Sidney Irrigated-Richland 0.0261606
## Sidney Irrigated-Sidney Dryland 0.0344455
nifH_irr_Observ<-ggboxplot(meta_nifH, x = "Plot", y = "Observed",
rug = TRUE,
fill = "Plot", xlab = " ",
palette = farm_col_paired)+
rremove("x.text")
nifH_irr_Shannon<-ggboxplot(meta_nifH, x = "Plot", y = "Shannon",
rug = TRUE,
fill = "Plot", xlab = " ",
palette = farm_col_paired)+
rremove("x.text")
nifH_irr_Chao<- ggboxplot(meta_nifH, x = "Plot", y = "Chao1",
rug = TRUE,
fill = "Plot", xlab = " ",
palette = farm_col_paired)+
rremove("x.text")
nifH_irr_InvSim<- ggboxplot(meta_nifH, x = "Plot", y = "Simpson",
rug = TRUE,
fill = "Plot", xlab = " ",
palette = farm_col_paired)+
rremove("x.text")
alpha_nifH_irr_fig<-ggarrange(nifH_irr_Observ, nifH_irr_Shannon, nifH_irr_Chao, nifH_irr_InvSim, ncol = 2, nrow = 2, common.legend = TRUE, legend = "bottom")
annotate_figure(alpha_nifH_irr_fig, top = text_grob("Alpha Diversity of nifH by Irrigation", size = 20))
aov_observed_irr_nifH <- aov(Observed ~ Plot, meta_nifH)
summary(aov_observed_irr_nifH)
## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 1351360 1351360 81.29 3.28e-12 ***
## Residuals 52 864472 16624
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
aov_observed_irr_nifH <- aov(Observed ~ Plot/Site, meta_nifH)
summary(aov_observed_irr_nifH)
## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 1351360 1351360 337.86 < 2e-16 ***
## Plot:Site 7 684482 97783 24.45 2.32e-13 ***
## Residuals 45 179990 4000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
aov_Chao1_irr_nifH <- aov(Chao1 ~ Plot, meta_nifH)
summary(aov_Chao1_irr_nifH)
## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 1328751 1328751 62.35 1.86e-10 ***
## Residuals 52 1108257 21313
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
aov_Chao1_irr_nifH <- aov(Chao1 ~ Plot/Site, meta_nifH)
summary(aov_Chao1_irr_nifH)
## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 1328751 1328751 234.21 < 2e-16 ***
## Plot:Site 7 852953 121850 21.48 2.11e-12 ***
## Residuals 45 255304 5673
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
aov_Simpson_irr_nifH <- aov(Simpson ~ Plot, meta_nifH)
summary(aov_Simpson_irr_nifH)
## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 0.1043 0.10431 19.64 4.85e-05 ***
## Residuals 52 0.2762 0.00531
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
aov_Simpson_irr_nifH <- aov(Simpson ~ Plot/Site, meta_nifH)
summary(aov_Simpson_irr_nifH)
## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 0.10431 0.10431 23.374 1.59e-05 ***
## Plot:Site 7 0.07539 0.01077 2.413 0.0347 *
## Residuals 45 0.20083 0.00446
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nifH_irr_Observ_scatter<-ggscatter(meta_nifH, x = "total_precip_irr", y = "Observed",
rug = TRUE,
add = "reg.line",
color = "Plot", xlab = " ",
palette = farm_col_paired)+
rremove("x.text")
nifH_irr_Shannon_scatter<-ggscatter(meta_nifH, x = "total_precip_irr", y = "Shannon",
rug = TRUE,
add = "reg.line",
color = "Plot", xlab = " ",
palette = farm_col_paired)+
rremove("x.text")
nifH_irr_Chao_scatter<- ggscatter(meta_nifH, x = "total_precip_irr", y = "Chao1",
rug = TRUE,
add = "reg.line",
color = "Plot", xlab = " ",
palette = farm_col_paired)+
rremove("x.text")
nifH_irr_InvSim_scatter<- ggscatter(meta_nifH, x = "total_precip_irr", y = "Simpson",
rug = TRUE,
add = "reg.line",
color = "Plot", xlab = " ",
palette = farm_col_paired)+
rremove("x.text")
alpha_nifH_irr_fig<-ggarrange(nifH_irr_Observ_scatter, nifH_irr_Shannon_scatter,
nifH_irr_Chao_scatter, nifH_irr_InvSim_scatter,
ncol = 2, nrow = 2, common.legend = TRUE, legend = "bottom")
annotate_figure(alpha_nifH_irr_fig,
top = text_grob("Alpha Diversity of nifH by Irrigation", size = 20))
https://joey711.github.io/phyloseq/plot_ordination-examples.html
phynifH_ord_RDA <- ordinate(physeq_nifH_ord, "RDA", "bray")
plot_ordination(physeq_nifH_ord, phynifH_ord_RDA, color = "Site", shape = "Plot")+
geom_point(size = 3)+
scale_color_manual(values = farm_col_paired)+
stat_ellipse(type = "norm", linetype = 2, aes(color ="Plot"), show.legend = TRUE) +
ggtitle("nifH RDA Ordination by Site and Irrigation")+
theme_bw()
Publish to the tiff
## png
## 2
plot_ordination(physeq_nifH_ord, phynifH_ord_RDA, color = "prev_crop", shape = "Plot")+
geom_point(size = 3)+
scale_color_manual(values = farm_col_paired)+
ggtitle("nifH RDA Ordination by Previous Crop and Irrigation")+
theme_bw()
plot_ordination(physeq_nifH_ord, phynifH_ord_RDA, color = "Tillage", shape = "Plot")+
geom_point(size = 3)+
scale_color_manual(values = farm_col_paired)+
ggtitle("nifH RDA Ordination by Previous Crop and Irrigation")+
theme_bw()
plot_ordination(physeq_nifH_ord, phynifH_ord_RDA, color = "Pea_variety", shape = "Plot")+
geom_point(size = 3)+
scale_color_manual(values = farm_col_paired)+
ggtitle("nifH RDA Ordination by Previous Crop and Irrigation")+
theme_bw()
####RDA has alot of overlap
phynifH_ord_NMDS <- ordinate(physeq_nifH_ord, "NMDS", "bray")
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1618779
## Run 1 stress 0.1666296
## Run 2 stress 0.1618753
## ... New best solution
## ... Procrustes: rmse 0.001579567 max resid 0.008369045
## ... Similar to previous best
## Run 3 stress 0.1887229
## Run 4 stress 0.161878
## ... Procrustes: rmse 0.001580814 max resid 0.008390491
## ... Similar to previous best
## Run 5 stress 0.1618754
## ... Procrustes: rmse 4.945684e-05 max resid 0.0003091343
## ... Similar to previous best
## Run 6 stress 0.1618781
## ... Procrustes: rmse 0.001584558 max resid 0.008371028
## ... Similar to previous best
## Run 7 stress 0.1657018
## Run 8 stress 0.1618754
## ... Procrustes: rmse 3.746306e-05 max resid 0.0002276462
## ... Similar to previous best
## Run 9 stress 0.1887328
## Run 10 stress 0.1618754
## ... Procrustes: rmse 9.792843e-05 max resid 0.0005519838
## ... Similar to previous best
## Run 11 stress 0.1618754
## ... Procrustes: rmse 7.118828e-05 max resid 0.0004423043
## ... Similar to previous best
## Run 12 stress 0.1618779
## ... Procrustes: rmse 0.001579864 max resid 0.008364055
## ... Similar to previous best
## Run 13 stress 0.1618779
## ... Procrustes: rmse 0.00158213 max resid 0.008386758
## ... Similar to previous best
## Run 14 stress 0.1618781
## ... Procrustes: rmse 0.001578848 max resid 0.008382921
## ... Similar to previous best
## Run 15 stress 0.1618779
## ... Procrustes: rmse 0.001582478 max resid 0.008396537
## ... Similar to previous best
## Run 16 stress 0.1618754
## ... Procrustes: rmse 4.310258e-05 max resid 0.0002645992
## ... Similar to previous best
## Run 17 stress 0.1618754
## ... Procrustes: rmse 9.174953e-05 max resid 0.0004160894
## ... Similar to previous best
## Run 18 stress 0.1618754
## ... Procrustes: rmse 4.384431e-05 max resid 0.0002702027
## ... Similar to previous best
## Run 19 stress 0.1618753
## ... Procrustes: rmse 2.118693e-05 max resid 0.00010052
## ... Similar to previous best
## Run 20 stress 0.1618754
## ... Procrustes: rmse 2.938269e-05 max resid 0.0001766266
## ... Similar to previous best
## *** Solution reached
plot_ordination(physeq_nifH_ord, phynifH_ord_NMDS)
phynifH_ord_NMDS
##
## Call:
## metaMDS(comm = veganifyOTU(physeq), distance = distance)
##
## global Multidimensional Scaling using monoMDS
##
## Data: wisconsin(sqrt(veganifyOTU(physeq)))
## Distance: bray
##
## Dimensions: 2
## Stress: 0.1618753
## Stress type 1, weak ties
## Two convergent solutions found after 20 tries
## Scaling: centring, PC rotation, halfchange scaling
## Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))'
####Plot NMDS with Site and Irrigation
plot_ordination(physeq_nifH_ord, phynifH_ord_NMDS, color = "Site")+
geom_point(size = 3)+
scale_color_manual(values = farm_col_paired)+
ggtitle("nifH NMDS Ordination by Site")+
theme_bw()
#export physeq OTU for network analysis and saving
nifH_trim_OTU <- as(otu_table(physeq_nifH_ord), "matrix")
ID <- rownames(nifH_trim_OTU)
nifH_trim_OTU_1 <- cbind(ID, nifH_trim_OTU)
rownames(nifH_trim_OTU_1)<-c()
write.table(nifH_trim_OTU_1, file = "nifH_trim_OTU.txt", sep = "\t", row.names = FALSE)
plot_ordination(physeq_nifH_ord, phynifH_ord_NMDS, color = "Pea_variety")+
geom_point(size = 3)+
scale_color_manual(values = farm_col_paired)+
ggtitle("nifH NMDS Ordination by Pea Variety")+
theme_bw()
plot_ordination(physeq_nifH_ord, phynifH_ord_NMDS, shape = "Plot", color = "Site")+
geom_point(size = 3)+
scale_color_manual(values = farm_col_paired)+
ggtitle("nifH NMDS Ordination by Irrigation method")+
stat_ellipse(type = "norm", linetype = 2, aes(color = "Plot")) +
theme_bw()
Publish to tiff
## png
## 2
plot_ordination(physeq_nifH_ord, phynifH_ord_NMDS, color = "total_precip_irr")+
geom_point(size = 3)+
scale_color_gradient(low='#05D9F6', high='#5011D1')+
ggtitle("nifH NMDS Ordination by Total Precipitation")+
theme_bw()
plot_ordination(physeq_nifH_ord, phynifH_ord_NMDS, color = "Moisture_Content")+
geom_point(size = 3)+
scale_color_continuous(low='#05D9F6', high='#5011D1')+
ggtitle("nifH NMDS Ordination by Mositure Content")+
theme_bw()
plot_ordination(physeq_nifH_ord, phynifH_ord_NMDS,color = "Tillage")+
geom_point(size = 3)+
scale_color_manual(values = farm_col_dark, breaks=c("Conventional", "Culti-roller", "No_till"),
labels=c("Conventional", "Culti Roller", "No Till"))+
ggtitle("nifH NMDS Ordination by Tillage")+
theme_bw()
Publish to .tiff
## png
## 2
plot_ordination(physeq_nifH_ord, phynifH_ord_NMDS, color = "prev_crop")+
geom_point(size = 3)+
scale_color_manual(values = c("#1B9E77", "#D95F02", "#7570B3", "#66a61E"),
name = "Previous Crop",
breaks=c("barley", "Chem_fallow", "Spring_wheat", "winter_wheat"),
labels=c("Barley", "Chemical Fallow", "Spring Wheat", "Winter Wheat"))+
ggtitle("nifH NMDS Ordination by Previous Crop")+
theme_bw()
Publish to .tiff
## png
## 2
####If a group (Site) in the MDS space are close but have diffrenent dispersion you could have a signifcant results when it is only a diffrence in dispersion.
Anderson (2006)-https://www.ncbi.nlm.nih.gov/pubmed/16542252 https://onlinelibrary.wiley.com/doi/epdf/10.1111/j.1461-0248.2006.00926.x
disp.plot <- betadisper(distance(physeq_nifH_ord, method = "bray"), meta2$Plot)
permutest(disp.plot, pairwise=TRUE, permutations=1000)
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 1000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.023867 0.0238672 6.7103 1000 0.01199 *
## Residuals 52 0.184953 0.0035568
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## dryland irrigated
## dryland 0.012
## irrigated 0.012408
boxplot(disp.plot)
plot(disp.plot, hull = FALSE, ellipse = TRUE)
disp.Till <- betadisper(distance(physeq_nifH_ord, method = "bray"), meta2$Tillage)
permutest(disp.Till, pairwise=TRUE, permutations=1000)
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 1000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.039609 0.0198044 4.2834 1000 0.01898 *
## Residuals 51 0.235799 0.0046235
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## Conventional Culti-roller No_till
## Conventional 9.9900e-04 0.0769
## Culti-roller 4.4826e-05 0.1399
## No_till 8.3245e-02 1.4099e-01
boxplot(disp.Till)
plot(disp.Till, hull = FALSE, ellipse = TRUE)
disp.prev_crop <- betadisper(distance(physeq_nifH_ord, method = "bray"), meta2$prev_crop)
permutest(disp.prev_crop, pairwise=TRUE, permutations=1000)
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 1000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.34683 0.115610 14.42 1000 0.000999 ***
## Residuals 50 0.40086 0.008017
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## barley Chem_fallow Spring_wheat winter_wheat
## barley 8.2917e-02 9.9900e-04 0.0010
## Chem_fallow 8.1169e-02 2.9970e-03 0.0020
## Spring_wheat 3.1442e-04 6.5930e-04 0.6963
## winter_wheat 2.7285e-05 3.6151e-05 6.7668e-01
boxplot(disp.prev_crop)
plot(disp.prev_crop, hull = FALSE, ellipse = TRUE)
TukeyHSD(disp.prev_crop)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = distances ~ group, data = df)
##
## $group
## diff lwr upr p adj
## Chem_fallow-barley -0.03352468 -0.1077209 0.04067151 0.6292204
## Spring_wheat-barley -0.18749276 -0.2996669 -0.07531866 0.0002814
## winter_wheat-barley -0.22719822 -0.3393723 -0.11502412 0.0000114
## Spring_wheat-Chem_fallow -0.15396808 -0.2625802 -0.04535597 0.0023889
## winter_wheat-Chem_fallow -0.19367354 -0.3022856 -0.08506143 0.0001047
## winter_wheat-Spring_wheat -0.03970546 -0.1770901 0.09767919 0.8684302
disp.site <- betadisper(distance(physeq_nifH_ord, method = "bray"), meta2$Site)
permutest(disp.site, pairwise=TRUE, permutations=1000)
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 1000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 8 0.35944 0.04493 3.2652 1000 0.005994 **
## Residuals 45 0.61921 0.01376
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## Conrad Corvallis Huntley Dryland Huntley Irrigated
## Conrad 0.03496503 0.03296703 0.06093906
## Corvallis 0.03646739 0.68131868 0.46953047
## Huntley Dryland 0.03003786 0.67097611 0.29470529
## Huntley Irrigated 0.06780670 0.44317758 0.29284235
## Kalispell 0.91870043 0.00041400 0.00027429 0.00246753
## Mocassin 0.52173883 0.07999231 0.06329790 0.16335615
## Richland 0.57388386 0.00439742 0.00283358 0.02313450
## Sidney Dryland 0.99250039 0.00143581 0.00099248 0.00635063
## Sidney Irrigated 0.32390060 0.21919454 0.17927660 0.38365390
## Kalispell Mocassin Richland Sidney Dryland
## Conrad 0.92007992 0.54045954 0.55844156 0.98601399
## Corvallis 0.00099900 0.08691309 0.00899101 0.00399600
## Huntley Dryland 0.00199800 0.07492507 0.00599401 0.00399600
## Huntley Irrigated 0.00899101 0.17082917 0.02897103 0.00699301
## Kalispell 0.31168831 0.26973027 0.86113886
## Mocassin 0.32083698 0.81818182 0.39160839
## Richland 0.27360130 0.82416350 0.37962038
## Sidney Dryland 0.86114730 0.40276122 0.39167738
## Sidney Irrigated 0.15341808 0.67667507 0.47279724 0.20368026
## Sidney Irrigated
## Conrad 0.3167
## Corvallis 0.2078
## Huntley Dryland 0.1738
## Huntley Irrigated 0.3956
## Kalispell 0.1409
## Mocassin 0.6643
## Richland 0.4446
## Sidney Dryland 0.2138
## Sidney Irrigated
boxplot(disp.site)
plot(disp.site, hull = FALSE, ellipse = TRUE)
TukeyHSD(disp.site)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = distances ~ group, data = df)
##
## $group
## diff lwr upr
## Corvallis-Conrad 0.1970866601 -0.02350369 0.417677013
## Huntley Dryland -Conrad 0.2056910662 -0.01489929 0.426281420
## Huntley Irrigated-Conrad 0.1730650455 -0.04752531 0.393655399
## Kalispell-Conrad -0.0092615904 -0.22985194 0.211328763
## Mocassin-Conrad 0.0683743608 -0.15221599 0.288964714
## Richland-Conrad 0.0514247120 -0.16916564 0.272015065
## Sidney Dryland-Conrad 0.0008764728 -0.21971388 0.221466826
## Sidney Irrigated-Conrad 0.1080798206 -0.11251053 0.328670174
## Huntley Dryland -Corvallis 0.0086044061 -0.21198595 0.229194759
## Huntley Irrigated-Corvallis -0.0240216145 -0.24461197 0.196568739
## Kalispell-Corvallis -0.2063482505 -0.42693860 0.014242103
## Mocassin-Corvallis -0.1287122993 -0.34930265 0.091878054
## Richland-Corvallis -0.1456619481 -0.36625230 0.074928405
## Sidney Dryland-Corvallis -0.1962101873 -0.41680054 0.024380166
## Sidney Irrigated-Corvallis -0.0890068395 -0.30959719 0.131583514
## Huntley Irrigated-Huntley Dryland -0.0326260207 -0.25321637 0.187964333
## Kalispell-Huntley Dryland -0.2149526566 -0.43554301 0.005637697
## Mocassin-Huntley Dryland -0.1373167054 -0.35790706 0.083273648
## Richland-Huntley Dryland -0.1542663542 -0.37485671 0.066323999
## Sidney Dryland-Huntley Dryland -0.2048145934 -0.42540495 0.015775760
## Sidney Irrigated-Huntley Dryland -0.0976112456 -0.31820160 0.122979108
## Kalispell-Huntley Irrigated -0.1823266359 -0.40291699 0.038263717
## Mocassin-Huntley Irrigated -0.1046906847 -0.32528104 0.115899669
## Richland-Huntley Irrigated -0.1216403336 -0.34223069 0.098950020
## Sidney Dryland-Huntley Irrigated -0.1721885728 -0.39277893 0.048401781
## Sidney Irrigated-Huntley Irrigated -0.0649852250 -0.28557558 0.155605128
## Mocassin-Kalispell 0.0776359512 -0.14295440 0.298226305
## Richland-Kalispell 0.0606863024 -0.15990405 0.281276656
## Sidney Dryland-Kalispell 0.0101380632 -0.21045229 0.230728417
## Sidney Irrigated-Kalispell 0.1173414110 -0.10324894 0.337931764
## Richland-Mocassin -0.0169496488 -0.23754000 0.203640704
## Sidney Dryland-Mocassin -0.0674978880 -0.28808824 0.153092465
## Sidney Irrigated-Mocassin 0.0397054598 -0.18088489 0.260295813
## Sidney Dryland-Richland -0.0505482392 -0.27113859 0.170042114
## Sidney Irrigated-Richland 0.0566551086 -0.16393524 0.277245462
## Sidney Irrigated-Sidney Dryland 0.1072033478 -0.11338701 0.327793701
## p adj
## Corvallis-Conrad 0.1134000
## Huntley Dryland -Conrad 0.0849905
## Huntley Irrigated-Conrad 0.2342545
## Kalispell-Conrad 1.0000000
## Mocassin-Conrad 0.9829927
## Richland-Conrad 0.9974322
## Sidney Dryland-Conrad 1.0000000
## Sidney Irrigated-Conrad 0.8019371
## Huntley Dryland -Corvallis 1.0000000
## Huntley Irrigated-Corvallis 0.9999912
## Kalispell-Corvallis 0.0830921
## Mocassin-Corvallis 0.6169037
## Richland-Corvallis 0.4537416
## Sidney Dryland-Corvallis 0.1166882
## Sidney Irrigated-Corvallis 0.9217819
## Huntley Irrigated-Huntley Dryland 0.9999076
## Kalispell-Huntley Dryland 0.0613855
## Mocassin-Huntley Dryland 0.5332525
## Richland-Huntley Dryland 0.3765001
## Sidney Dryland-Huntley Dryland 0.0875794
## Sidney Irrigated-Huntley Dryland 0.8754541
## Kalispell-Huntley Irrigated 0.1797308
## Mocassin-Huntley Irrigated 0.8276676
## Richland-Huntley Irrigated 0.6843270
## Sidney Dryland-Huntley Irrigated 0.2399571
## Sidney Irrigated-Huntley Irrigated 0.9876732
## Mocassin-Kalispell 0.9633505
## Richland-Kalispell 0.9921001
## Sidney Dryland-Kalispell 1.0000000
## Sidney Irrigated-Kalispell 0.7237181
## Richland-Mocassin 0.9999994
## Sidney Dryland-Mocassin 0.9843152
## Sidney Irrigated-Mocassin 0.9995991
## Sidney Dryland-Richland 0.9977226
## Sidney Irrigated-Richland 0.9950077
## Sidney Irrigated-Sidney Dryland 0.8087566
####The results from the betadispersion show that we have a signficant diffrence in the heterogeneity of our sites due to each of the farm amnagament factors. We see that there are sites that do not have signifcant diffrences and others that do. Beta dispersion cannot except models so we cannot detect if the diffrences are due to nestedness (prev_crop/Site).
####Do to the design of the experiment it will be hard to determine if the farm managment practices are responsible for the variation in bacterial population.
####We will test the significance of the farm managment using the permutated-mulitvariate-ANOVA function in vegan called adonis. Adonis can test models and nestedness though the above beta dispersion test showed that most of the signficance is due to the diseprsion we cannot say much. But we can show we have the location effect.
adonis(distance(physeq_nifH_ord, method = "bray") ~Plot*prev_crop*Tillage*Pea_variety, data = meta2, permutations = 1000)
##
## Call:
## adonis(formula = distance(physeq_nifH_ord, method = "bray") ~ Plot * prev_crop * Tillage * Pea_variety, data = meta2, permutations = 1000)
##
## Permutation: free
## Number of permutations: 1000
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Plot 1 1.4006 1.40059 3.8638 0.07104 0.000999 ***
## prev_crop 2 2.5341 1.26705 3.4954 0.12853 0.000999 ***
## Tillage 2 3.0720 1.53601 4.2374 0.15581 0.000999 ***
## Pea_variety 6 1.3480 0.22466 0.6198 0.06837 1.000000
## Plot:Tillage 1 1.5247 1.52471 4.2062 0.07733 0.000999 ***
## Plot:Pea_variety 6 0.9169 0.15282 0.4216 0.04651 1.000000
## prev_crop:Pea_variety 9 1.7119 0.19021 0.5247 0.08683 1.000000
## Tillage:Pea_variety 10 2.1740 0.21740 0.5997 0.11026 1.000000
## Plot:Tillage:Pea_variety 4 0.6841 0.17103 0.4718 0.03470 0.999001
## Residuals 12 4.3499 0.36249 0.22063
## Total 53 19.7163 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis(distance(physeq_nifH_ord, method = "bray") ~(Plot*prev_crop*Tillage*Pea_variety)/Site, data = meta2, permutations = 1000)
##
## Call:
## adonis(formula = distance(physeq_nifH_ord, method = "bray") ~ (Plot * prev_crop * Tillage * Pea_variety)/Site, data = meta2, permutations = 1000)
##
## Permutation: free
## Number of permutations: 1000
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model
## Plot 1 1.4006 1 0
## prev_crop 2 2.5341 1 0
## Tillage 2 3.0720 2 0
## Pea_variety 6 1.3480 0 0
## Plot:Tillage 1 1.5247 2 0
## Plot:Pea_variety 6 0.9169 0 0
## prev_crop:Pea_variety 9 1.7119 0 0
## Tillage:Pea_variety 10 2.1740 0 0
## Plot:Tillage:Pea_variety 4 0.6841 0 0
## Plot:prev_crop:Tillage:Pea_variety:Site 12 4.3499 0 0
## Residuals 0 0.0000 -Inf
## Total 53 19.7163
## R2 Pr(>F)
## Plot 0.07104 1
## prev_crop 0.12853 1
## Tillage 0.15581 1
## Pea_variety 0.06837 1
## Plot:Tillage 0.07733 1
## Plot:Pea_variety 0.04651 1
## prev_crop:Pea_variety 0.08683 1
## Tillage:Pea_variety 0.11026 1
## Plot:Tillage:Pea_variety 0.03470 1
## Plot:prev_crop:Tillage:Pea_variety:Site 0.22063 1
## Residuals 0.00000
## Total 1.00000
####Since there is issues with doing permutated anova over multivatrite data lets try to fit the chemical and farm managment data the NMDS orientation space using the envfit function in vegan
Model Selection
meta3<-meta2[,-c(3,4,10,11,27,29,35,38,42,46)]
meta3<-meta3[,-c(2,11:13)]
####Removing Sulfate_Sulfur, Boron, Molybdenum, Potassium, Vanadium, Chromium and Sodium (From chemical_analysis.rmd)
meta3<-meta3[,-c(16,18,21,29,31)]
envfitnifH <- envfit(phynifH_ord_NMDS , meta3, na.rm = TRUE, permu= 999)
envfitnifH
##
## ***VECTORS
##
## NMDS1 NMDS2 r2 Pr(>r)
## season_precip 0.30096 -0.95364 0.1693 0.008 **
## irrgation 0.54023 0.84152 0.1767 0.007 **
## total_precip_irr 0.79985 -0.60021 0.1800 0.008 **
## grain_yield -0.10337 0.99464 0.0052 0.870
## Organic_Matter -0.58271 -0.81268 0.5434 0.001 ***
## Moisture_Content 0.21634 -0.97632 0.2847 0.002 **
## Nitrate_Nitrite 0.01586 0.99987 0.0848 0.092 .
## Ammonia -0.16640 -0.98606 0.0684 0.163
## Av_Phosphorus 0.77470 0.63232 0.2592 0.001 ***
## Av_Potassium 0.98960 -0.14383 0.3455 0.001 ***
## pH -0.02160 -0.99977 0.0775 0.123
## Barium -0.09002 -0.99594 0.3593 0.001 ***
## Calcium 0.29225 -0.95634 0.2220 0.003 **
## Cobalt 0.00453 -0.99999 0.0831 0.115
## Copper 0.31097 -0.95042 0.1409 0.022 *
## Iron 0.25017 -0.96820 0.1241 0.031 *
## Magnesium 0.22694 -0.97391 0.0394 0.342
## Manganese 0.18760 -0.98225 0.1245 0.045 *
## Nickel 0.63393 -0.77339 0.1124 0.033 *
## Phosphorus 0.09644 -0.99534 0.2283 0.004 **
## Sulfur 0.10985 -0.99395 0.3851 0.001 ***
## Zinc 0.44396 -0.89605 0.2825 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
##
## ***FACTORS:
##
## Centroids:
## NMDS1 NMDS2
## SiteConrad 0.9037 0.2518
## SiteCorvallis 0.0861 0.7313
## SiteHuntley Dryland -0.2456 0.0338
## SiteHuntley Irrigated -0.1426 -0.0653
## SiteKalispell -1.6148 -0.2105
## SiteMocassin 0.0258 -0.7600
## SiteRichland 0.2990 -0.0971
## SiteSidney Dryland -0.1793 0.4292
## SiteSidney Irrigated 0.8678 -0.3131
## Pea_varietyAC Earlystar -0.1321 0.1464
## Pea_varietyCDC Saffron -0.0296 -0.0669
## Pea_varietyCDC Saffron 0.4876 0.1231
## Pea_varietyDelta 0.1280 -0.0638
## Pea_varietyDS Admiral 0.0304 -0.0922
## Pea_varietyMajoret -0.0978 0.0629
## Pea_varietyNavarro -0.0138 -0.0287
## Plotdryland 0.1607 -0.0285
## Plotirrigated -0.2009 0.0356
## TillageConventional -0.3088 -0.0315
## TillageCulti-roller 0.0861 0.7313
## TillageNo_till 0.1681 -0.1274
## prev_cropbarley -0.5571 0.1518
## prev_cropChem_fallow 0.1944 0.1544
## prev_cropSpring_wheat 0.8678 -0.3131
## prev_cropwinter_wheat 0.0258 -0.7600
##
## Goodness of fit:
## r2 Pr(>r)
## Site 0.8843 0.001 ***
## Pea_variety 0.0323 0.991
## Plot 0.0449 0.107
## Tillage 0.1578 0.005 **
## prev_crop 0.4008 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
plot(phynifH_ord_NMDS, display = "sites")
plot(envfitnifH, p.max = 0.001 )
“The function calculates a community dissimilarity matrix using vegdist. Then it selects all possible subsets of environmental variables, scales the variables, and calculates Euclidean distances for this subset using dist. The function finds the correlation between community dissimilarities and environmental distances, and for each size of subsets, saves the best result. There are 2^p-1 subsets of p variables, and an exhaustive search may take a very, very, very long time (parameter upto offers a partial relief).”
“Mantel statistic is simply a correlation between entries of two dissimilarity matrices (some use cross products, but these are linearly related). However, the significance cannot be directly assessed, because there are N(N-1)/2 entries for just N observations. Mantel developed asymptotic test, but here we use permutations of N rows and columns of dissimilarity matrix.”
OTU_nifH_trim<- as(otu_table(physeq_nifH_ord ), "matrix")
#Transpose the data to have sample names on rows
abund_tablenifH<-t(OTU_nifH_trim)
nrow(abund_tablenifH)
## [1] 54
setdiff(rownames(meta3), rownames(abund_tablenifH))
## character(0)
#First detect amount of cores avalible
detectCores()
## [1] 12
#get bray-curtis distance
abund_distnifH<-vegdist(abund_tablenifH, method = "bray")
#make bioenv model against whole meta table with ".", and use gower metric to measure distance in order to incroprate the factor variables. Model will be used up to 7 vriables in final model. parelle processing on 8 cores
#nifH.bioenv <- bioenv(abund_distnifH ~ ., meta3, index = "bray", method = "pearson",
# metric = "gower", upto = 10, parallel = 10)
#summary(nifH.bioenv)
#nifH.bioenvdist<-bioenvdist(nifH.bioenv, which = "best")
#mantel(nifH.bioenvdist, abund_distnifH)
#nifH.bioenv$whichbest
####The best model from the inital BioEnv model selection shows that the Site explains 59% of the variation in the bray-curitis distance. But the best model explains 63% of the variation when Organic Matter and nitritie and nitrate are added in.
meta4<-meta3[,-c(1)]
#nifH.bioenv.site.na <- bioenv(abund_distnifH ~ ., meta4, index = "bray",
# method = "pearson", metric = "gower", upto = 10, parallel = 10)
#summary(nifH.bioenv.site.na)
#nifH.bioenvdist.site.na<-bioenvdist(nifH.bioenv.site.na, which = "best")
#mantel(nifH.bioenvdist.site.na, abund_distnifH)
#nifH.bioenv.site.na$whichbest
#nifH.bioenv.chem <- bioenv(abund_tablenifH ~ Organic_Matter + Moisture_Content +
# Nitrate_Nitrite + Ammonia + Av_Phosphorus + Av_Potassium +
# pH + Barium + Calcium + Cobalt + Copper + Iron + Magnesium +
# Manganese + Nickel + Phosphorus + Sulfur + Zinc, meta3,
# index = "bray", method = "pearson", metric = "gower",
# upto = 10, parallel = 10)
#summary(nifH.bioenv.chem)
#nifH.bioenvdist.chem<-bioenvdist(nifH.bioenv.site.na, which = "best")
#mantel(nifH.bioenvdist.site.na, abund_distnifH)
#nifH.bioenv.chem$whichbest
m1_nifH <- cca(abund_tablenifH ~ ., meta3)
m0_nifH <- cca(abund_tablenifH ~ 1, meta3)
m1_nifH
## Call: cca(formula = abund_tablenifH ~ Site + Pea_variety + Plot +
## season_precip + irrgation + total_precip_irr + Tillage + prev_crop
## + grain_yield + Organic_Matter + Moisture_Content +
## Nitrate_Nitrite + Ammonia + Av_Phosphorus + Av_Potassium + pH +
## Barium + Calcium + Cobalt + Copper + Iron + Magnesium + Manganese
## + Nickel + Phosphorus + Sulfur + Zinc, data = meta3)
##
## Inertia Proportion Rank
## Total 10.7659 1.0000
## Constrained 7.8889 0.7328 33
## Unconstrained 2.8771 0.2672 20
## Inertia is scaled Chi-square
## Some constraints were aliased because they were collinear (redundant)
##
## Eigenvalues for constrained axes:
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6 CCA7 CCA8 CCA9 CCA10
## 0.8070 0.6757 0.6117 0.5408 0.5291 0.4814 0.4024 0.3842 0.3332 0.3174
## CCA11 CCA12 CCA13 CCA14 CCA15 CCA16 CCA17 CCA18 CCA19 CCA20
## 0.2820 0.2391 0.2262 0.2006 0.1894 0.1707 0.1638 0.1537 0.1350 0.1302
## CCA21 CCA22 CCA23 CCA24 CCA25 CCA26 CCA27 CCA28 CCA29 CCA30
## 0.1147 0.1077 0.0963 0.0872 0.0804 0.0770 0.0695 0.0613 0.0515 0.0497
## CCA31 CCA32 CCA33
## 0.0441 0.0429 0.0330
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.29098 0.28225 0.24918 0.23022 0.20041 0.17909 0.17644 0.16025
## (Showed only 8 of all 20 unconstrained eigenvalues)
m0_nifH
## Call: cca(formula = abund_tablenifH ~ 1, data = meta3)
##
## Inertia Rank
## Total 10.77
## Unconstrained 10.77 53
## Inertia is scaled Chi-square
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.8189 0.6921 0.6401 0.5860 0.5584 0.5037 0.4818 0.4700
## (Showed only 8 of all 53 unconstrained eigenvalues)
model_nifH <-ordistep(m0_nifH, scope=formula(m1_nifH))
##
## Start: abund_tablenifH ~ 1
##
## Df AIC F Pr(>F)
## + Site 8 868.86 3.0678 0.005 **
## + Tillage 2 873.90 3.2418 0.005 **
## + Organic_Matter 1 874.75 3.5957 0.005 **
## + prev_crop 3 875.14 2.3857 0.005 **
## + Nitrate_Nitrite 1 875.40 2.9283 0.005 **
## + irrgation 1 875.41 2.9174 0.005 **
## + season_precip 1 875.67 2.6597 0.005 **
## + Zinc 1 875.68 2.6463 0.005 **
## + Av_Potassium 1 875.71 2.6178 0.005 **
## + total_precip_irr 1 875.74 2.5881 0.005 **
## + Nickel 1 875.88 2.4512 0.005 **
## + Cobalt 1 875.92 2.4074 0.005 **
## + Barium 1 875.93 2.3929 0.005 **
## + Iron 1 875.93 2.3923 0.005 **
## + Copper 1 875.99 2.3375 0.005 **
## + Plot 1 876.00 2.3303 0.005 **
## + Magnesium 1 876.00 2.3264 0.005 **
## + Av_Phosphorus 1 876.01 2.3202 0.005 **
## + Phosphorus 1 876.07 2.2549 0.005 **
## + Moisture_Content 1 876.18 2.1464 0.005 **
## + Manganese 1 876.18 2.1462 0.005 **
## + grain_yield 1 876.27 2.0533 0.005 **
## + Ammonia 1 876.47 1.8528 0.005 **
## + Sulfur 1 876.47 1.8514 0.005 **
## + Calcium 1 876.57 1.7584 0.005 **
## + pH 1 876.84 1.4897 0.015 *
## + Pea_variety 6 882.82 0.8462 0.980
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Site
##
## Df AIC F Pr(>F)
## - Site 8 876.36 3.0678 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Nitrate_Nitrite 1 868.08 2.3207 0.005 **
## + Ammonia 1 869.09 1.4610 0.030 *
## + Av_Phosphorus 1 869.18 1.3900 0.060 .
## + Pea_variety 6 872.61 1.0730 0.185
## + Cobalt 1 869.64 1.0051 0.370
## + Manganese 1 869.65 0.9978 0.420
## + grain_yield 1 869.68 0.9722 0.445
## + Iron 1 869.74 0.9187 0.460
## + Phosphorus 1 869.79 0.8751 0.565
## + Copper 1 869.81 0.8624 0.565
## + Zinc 1 869.77 0.8984 0.575
## + Barium 1 869.78 0.8907 0.600
## + Av_Potassium 1 869.86 0.8240 0.705
## + Magnesium 1 869.87 0.8137 0.705
## + Nickel 1 869.82 0.8508 0.730
## + pH 1 869.81 0.8615 0.750
## + Organic_Matter 1 870.05 0.6648 0.900
## + Moisture_Content 1 869.98 0.7224 0.910
## + Calcium 1 870.08 0.6382 0.925
## + Sulfur 1 870.03 0.6823 0.950
## + Plot 0 868.86
## + season_precip 0 868.86
## + irrgation 0 868.86
## + total_precip_irr 0 868.86
## + Tillage 0 868.86
## + prev_crop 0 868.86
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Site + Nitrate_Nitrite
##
## Df AIC F Pr(>F)
## - Nitrate_Nitrite 1 868.86 2.3207 0.005 **
## - Site 8 875.40 2.9709 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Ammonia 1 868.19 1.5305 0.025 *
## + Av_Phosphorus 1 868.30 1.4419 0.070 .
## + Pea_variety 6 871.61 1.0759 0.190
## + Cobalt 1 868.83 1.0065 0.445
## + Zinc 1 868.92 0.9313 0.450
## + Manganese 1 868.81 1.0230 0.455
## + Iron 1 868.90 0.9530 0.495
## + Phosphorus 1 868.90 0.9475 0.525
## + Copper 1 868.95 0.9077 0.550
## + Barium 1 868.93 0.9309 0.580
## + Av_Potassium 1 868.97 0.8946 0.600
## + Magnesium 1 869.03 0.8495 0.665
## + pH 1 868.99 0.8775 0.690
## + Nickel 1 869.05 0.8322 0.765
## + grain_yield 1 869.10 0.7851 0.840
## + Organic_Matter 1 869.20 0.7088 0.875
## + Moisture_Content 1 869.15 0.7450 0.900
## + Sulfur 1 869.21 0.6989 0.920
## + Calcium 1 869.26 0.6602 0.925
## + Plot 0 868.08
## + season_precip 0 868.08
## + irrgation 0 868.08
## + total_precip_irr 0 868.08
## + Tillage 0 868.08
## + prev_crop 0 868.08
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Site + Nitrate_Nitrite + Ammonia
##
## Df AIC F Pr(>F)
## - Ammonia 1 868.08 1.5305 0.035 *
## - Nitrate_Nitrite 1 869.09 2.3726 0.005 **
## - Site 8 875.48 2.8982 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Av_Phosphorus 1 868.52 1.3192 0.080 .
## + Pea_variety 6 871.49 1.0780 0.155
## + Iron 1 868.97 0.9609 0.425
## + Cobalt 1 868.93 0.9935 0.470
## + Manganese 1 868.94 0.9858 0.475
## + Zinc 1 869.01 0.9343 0.505
## + Barium 1 869.00 0.9394 0.515
## + Copper 1 869.03 0.9155 0.535
## + Phosphorus 1 868.98 0.9541 0.545
## + pH 1 869.03 0.9173 0.620
## + Magnesium 1 869.10 0.8613 0.670
## + Av_Potassium 1 869.11 0.8515 0.720
## + Nickel 1 869.12 0.8419 0.740
## + grain_yield 1 869.16 0.8129 0.825
## + Organic_Matter 1 869.28 0.7172 0.885
## + Sulfur 1 869.27 0.7237 0.905
## + Calcium 1 869.33 0.6769 0.935
## + Moisture_Content 1 869.30 0.6981 0.960
## + Plot 0 868.19
## + season_precip 0 868.19
## + irrgation 0 868.19
## + total_precip_irr 0 868.19
## + Tillage 0 868.19
## + prev_crop 0 868.19
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1_nifH_site_na <- cca(abund_tablenifH ~ ., meta4)
m0_nifH_site_na <- cca(abund_tablenifH ~ 1, meta4)
m1_nifH_site_na
## Call: cca(formula = abund_tablenifH ~ Pea_variety + Plot +
## season_precip + irrgation + total_precip_irr + Tillage + prev_crop
## + grain_yield + Organic_Matter + Moisture_Content +
## Nitrate_Nitrite + Ammonia + Av_Phosphorus + Av_Potassium + pH +
## Barium + Calcium + Cobalt + Copper + Iron + Magnesium + Manganese
## + Nickel + Phosphorus + Sulfur + Zinc, data = meta4)
##
## Inertia Proportion Rank
## Total 10.7659 1.0000
## Constrained 7.6960 0.7148 32
## Unconstrained 3.0700 0.2852 21
## Inertia is scaled Chi-square
## Some constraints were aliased because they were collinear (redundant)
##
## Eigenvalues for constrained axes:
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6 CCA7 CCA8 CCA9 CCA10
## 0.8067 0.6725 0.6096 0.5394 0.5055 0.4714 0.3948 0.3835 0.3273 0.3114
## CCA11 CCA12 CCA13 CCA14 CCA15 CCA16 CCA17 CCA18 CCA19 CCA20
## 0.2705 0.2277 0.2256 0.1944 0.1815 0.1687 0.1585 0.1494 0.1316 0.1247
## CCA21 CCA22 CCA23 CCA24 CCA25 CCA26 CCA27 CCA28 CCA29 CCA30
## 0.1138 0.1030 0.0960 0.0871 0.0803 0.0698 0.0633 0.0572 0.0498 0.0445
## CCA31 CCA32
## 0.0430 0.0334
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.29949 0.28234 0.26627 0.23089 0.20393 0.19999 0.17647 0.17203
## (Showed only 8 of all 21 unconstrained eigenvalues)
m0_nifH_site_na
## Call: cca(formula = abund_tablenifH ~ 1, data = meta4)
##
## Inertia Rank
## Total 10.77
## Unconstrained 10.77 53
## Inertia is scaled Chi-square
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.8189 0.6921 0.6401 0.5860 0.5584 0.5037 0.4818 0.4700
## (Showed only 8 of all 53 unconstrained eigenvalues)
model_nifH <-ordistep(m0_nifH_site_na, scope=formula(m1_nifH_site_na))
##
## Start: abund_tablenifH ~ 1
##
## Df AIC F Pr(>F)
## + Tillage 2 873.90 3.2418 0.005 **
## + Organic_Matter 1 874.75 3.5957 0.005 **
## + prev_crop 3 875.14 2.3857 0.005 **
## + Nitrate_Nitrite 1 875.40 2.9283 0.005 **
## + irrgation 1 875.41 2.9174 0.005 **
## + season_precip 1 875.67 2.6597 0.005 **
## + Zinc 1 875.68 2.6463 0.005 **
## + Av_Potassium 1 875.71 2.6178 0.005 **
## + total_precip_irr 1 875.74 2.5881 0.005 **
## + Nickel 1 875.88 2.4512 0.005 **
## + Cobalt 1 875.92 2.4074 0.005 **
## + Barium 1 875.93 2.3929 0.005 **
## + Iron 1 875.93 2.3923 0.005 **
## + Copper 1 875.99 2.3375 0.005 **
## + Plot 1 876.00 2.3303 0.005 **
## + Magnesium 1 876.00 2.3264 0.005 **
## + Av_Phosphorus 1 876.01 2.3202 0.005 **
## + Phosphorus 1 876.07 2.2549 0.005 **
## + Moisture_Content 1 876.18 2.1464 0.005 **
## + Manganese 1 876.18 2.1462 0.005 **
## + grain_yield 1 876.27 2.0533 0.005 **
## + Ammonia 1 876.47 1.8528 0.005 **
## + Calcium 1 876.57 1.7584 0.005 **
## + Sulfur 1 876.47 1.8514 0.010 **
## + pH 1 876.84 1.4897 0.020 *
## + Pea_variety 6 882.82 0.8462 0.980
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Tillage
##
## Df AIC F Pr(>F)
## - Tillage 2 876.36 3.2418 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + prev_crop 3 871.82 2.5825 0.005 **
## + Organic_Matter 1 872.19 3.5560 0.005 **
## + total_precip_irr 1 872.60 3.1476 0.005 **
## + irrgation 1 872.87 2.8842 0.005 **
## + Av_Potassium 1 872.95 2.8096 0.005 **
## + Av_Phosphorus 1 872.98 2.7771 0.005 **
## + season_precip 1 873.16 2.6012 0.005 **
## + Magnesium 1 873.34 2.4237 0.005 **
## + Phosphorus 1 873.44 2.3340 0.005 **
## + Barium 1 873.62 2.1603 0.005 **
## + Nickel 1 873.63 2.1424 0.005 **
## + Moisture_Content 1 873.65 2.1231 0.005 **
## + Plot 1 873.71 2.0722 0.005 **
## + Manganese 1 873.73 2.0526 0.005 **
## + Nitrate_Nitrite 1 873.77 2.0084 0.005 **
## + Zinc 1 873.94 1.8456 0.005 **
## + Ammonia 1 873.95 1.8361 0.005 **
## + Calcium 1 874.00 1.7954 0.005 **
## + Iron 1 874.06 1.7324 0.005 **
## + Sulfur 1 874.08 1.7112 0.005 **
## + Copper 1 874.13 1.6675 0.005 **
## + grain_yield 1 874.17 1.6277 0.005 **
## + Cobalt 1 873.91 1.8798 0.015 *
## + pH 1 874.58 1.2350 0.130
## + Pea_variety 6 879.71 0.9109 0.850
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Tillage + prev_crop
##
## Df AIC F Pr(>F)
## - prev_crop 3 873.90 2.5825 0.005 **
## - Tillage 2 875.14 3.4833 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + irrgation 1 870.32 3.1471 0.005 **
## + Organic_Matter 1 870.67 2.8207 0.005 **
## + Ammonia 1 870.86 2.6513 0.005 **
## + total_precip_irr 1 871.34 2.2116 0.005 **
## + grain_yield 1 871.34 2.2108 0.005 **
## + Magnesium 1 871.34 2.2086 0.005 **
## + Barium 1 871.51 2.0507 0.005 **
## + season_precip 1 871.56 2.0064 0.005 **
## + Nickel 1 871.58 1.9892 0.005 **
## + Manganese 1 871.63 1.9468 0.005 **
## + pH 1 871.84 1.7541 0.005 **
## + Nitrate_Nitrite 1 871.63 1.9486 0.010 **
## + Cobalt 1 871.63 1.9418 0.010 **
## + Av_Potassium 1 871.94 1.6614 0.010 **
## + Calcium 1 872.01 1.6001 0.010 **
## + Zinc 1 871.95 1.6604 0.015 *
## + Copper 1 872.09 1.5325 0.020 *
## + Iron 1 871.81 1.7800 0.030 *
## + Av_Phosphorus 1 872.28 1.3598 0.045 *
## + Phosphorus 1 872.28 1.3633 0.120
## + Moisture_Content 1 872.41 1.2435 0.155
## + Sulfur 1 872.55 1.1177 0.340
## + Pea_variety 6 876.76 0.9773 0.570
## + Plot 0 871.82
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Tillage + prev_crop + irrgation
##
## Df AIC F Pr(>F)
## - irrgation 1 871.82 3.1471 0.005 **
## - prev_crop 3 872.87 2.6883 0.005 **
## - Tillage 2 872.96 3.0756 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Ammonia 1 870.03 1.9949 0.005 **
## + Manganese 1 870.12 1.9138 0.005 **
## + season_precip 1 870.13 1.9074 0.005 **
## + total_precip_irr 1 870.13 1.9074 0.005 **
## + grain_yield 1 870.20 1.8412 0.005 **
## + pH 1 870.27 1.7801 0.005 **
## + Av_Potassium 1 870.41 1.6562 0.005 **
## + Zinc 1 870.45 1.6182 0.005 **
## + Nitrate_Nitrite 1 870.01 2.0078 0.010 **
## + Organic_Matter 1 870.41 1.6546 0.010 **
## + Nickel 1 870.57 1.5191 0.010 **
## + Magnesium 1 870.54 1.5390 0.030 *
## + Copper 1 870.59 1.5016 0.040 *
## + Av_Phosphorus 1 870.70 1.3994 0.060 .
## + Cobalt 1 870.40 1.6632 0.090 .
## + Iron 1 870.65 1.4446 0.105
## + Moisture_Content 1 870.76 1.3526 0.115
## + Barium 1 870.95 1.1863 0.165
## + Phosphorus 1 870.83 1.2836 0.220
## + Sulfur 1 871.08 1.0655 0.375
## + Pea_variety 6 874.92 1.0036 0.520
## + Calcium 1 871.27 0.9011 0.685
## + Plot 0 870.32
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Tillage + prev_crop + irrgation + Ammonia
##
## Df AIC F Pr(>F)
## - Ammonia 1 870.32 1.9949 0.005 **
## - irrgation 1 870.86 2.4742 0.005 **
## - Tillage 2 872.15 2.7622 0.005 **
## - prev_crop 3 872.66 2.6590 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + season_precip 1 869.61 2.0562 0.005 **
## + total_precip_irr 1 869.61 2.0562 0.005 **
## + Nitrate_Nitrite 1 869.62 2.0545 0.005 **
## + pH 1 869.82 1.8745 0.005 **
## + Manganese 1 869.88 1.8270 0.005 **
## + Av_Potassium 1 870.01 1.7141 0.010 **
## + Organic_Matter 1 870.04 1.6910 0.010 **
## + Zinc 1 870.10 1.6350 0.010 **
## + grain_yield 1 870.11 1.6301 0.010 **
## + Nickel 1 870.13 1.6135 0.020 *
## + Copper 1 870.20 1.5505 0.025 *
## + Magnesium 1 870.33 1.4337 0.055 .
## + Iron 1 870.39 1.3872 0.135
## + Cobalt 1 870.43 1.3555 0.220
## + Av_Phosphorus 1 870.75 1.0745 0.250
## + Pea_variety 6 874.34 1.0205 0.365
## + Phosphorus 1 870.85 0.9908 0.435
## + Barium 1 870.85 0.9879 0.510
## + Sulfur 1 870.95 0.9063 0.625
## + Calcium 1 870.93 0.9229 0.655
## + Moisture_Content 1 871.02 0.8480 0.730
## + Plot 0 870.03
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Tillage + prev_crop + irrgation + Ammonia + season_precip
##
## Df AIC F Pr(>F)
## - season_precip 1 870.03 2.0562 0.005 **
## - Ammonia 1 870.13 2.1422 0.005 **
## - irrgation 1 870.58 2.5391 0.005 **
## - Tillage 2 871.49 2.5867 0.005 **
## - prev_crop 3 871.65 2.4051 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + grain_yield 1 869.54 1.7222 0.010 **
## + Nitrate_Nitrite 1 868.86 2.3031 0.020 *
## + Cobalt 1 869.93 1.3919 0.160
## + Manganese 1 869.86 1.4507 0.170
## + Av_Phosphorus 1 870.26 1.1212 0.270
## + Pea_variety 6 873.52 1.0515 0.310
## + Magnesium 1 870.30 1.0830 0.355
## + Iron 1 870.31 1.0795 0.360
## + pH 1 870.37 1.0279 0.415
## + Av_Potassium 1 870.39 1.0079 0.460
## + Phosphorus 1 870.39 1.0080 0.470
## + Barium 1 870.41 0.9929 0.515
## + Zinc 1 870.53 0.8910 0.575
## + Organic_Matter 1 870.50 0.9177 0.610
## + Nickel 1 870.60 0.8379 0.645
## + Sulfur 1 870.54 0.8863 0.660
## + Copper 1 870.61 0.8285 0.680
## + Moisture_Content 1 870.68 0.7700 0.800
## + Calcium 1 870.82 0.6544 0.915
## + Plot 0 869.61
## + total_precip_irr 0 869.61
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Tillage + prev_crop + irrgation + Ammonia + season_precip + grain_yield
##
## Df AIC F Pr(>F)
## - Ammonia 1 869.37 1.5161 0.015 *
## - grain_yield 1 869.61 1.7222 0.005 **
## - season_precip 1 870.11 2.1400 0.005 **
## - irrgation 1 870.60 2.5682 0.005 **
## - Tillage 2 871.60 2.6131 0.005 **
## - prev_crop 3 871.71 2.3952 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Nitrate_Nitrite 1 868.62 2.3907 0.010 **
## + Av_Phosphorus 1 870.05 1.2029 0.160
## + Pea_variety 6 873.33 1.0408 0.290
## + Cobalt 1 870.30 0.9985 0.415
## + Av_Potassium 1 870.34 0.9699 0.420
## + Iron 1 870.36 0.9536 0.425
## + Manganese 1 870.25 1.0415 0.465
## + pH 1 870.30 1.0004 0.505
## + Phosphorus 1 870.35 0.9605 0.535
## + Zinc 1 870.44 0.8877 0.540
## + Magnesium 1 870.46 0.8682 0.605
## + Barium 1 870.42 0.9015 0.620
## + Copper 1 870.49 0.8450 0.685
## + Nickel 1 870.49 0.8494 0.740
## + Sulfur 1 870.62 0.7421 0.830
## + Organic_Matter 1 870.62 0.7391 0.850
## + Moisture_Content 1 870.65 0.7188 0.890
## + Calcium 1 870.70 0.6781 0.905
## + Plot 0 869.54
## + total_precip_irr 0 869.54
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Tillage + prev_crop + irrgation + Ammonia + season_precip + grain_yield + Nitrate_Nitrite
##
## Df AIC F Pr(>F)
## - Ammonia 1 868.61 1.6173 0.015 *
## - grain_yield 1 868.86 1.8213 0.005 **
## - season_precip 1 869.33 2.2157 0.005 **
## - Nitrate_Nitrite 1 869.54 2.3907 0.005 **
## - irrgation 1 869.85 2.6529 0.005 **
## - Tillage 2 870.93 2.6640 0.005 **
## - prev_crop 3 871.16 2.4551 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Av_Phosphorus 1 869.14 1.1639 0.225
## + Pea_variety 6 872.23 1.0371 0.365
## + Cobalt 1 869.31 1.0309 0.370
## + Manganese 1 869.27 1.0642 0.380
## + Av_Potassium 1 869.31 1.0299 0.390
## + pH 1 869.33 1.0128 0.435
## + Phosphorus 1 869.31 1.0293 0.450
## + Iron 1 869.37 0.9871 0.485
## + Zinc 1 869.47 0.9046 0.525
## + Magnesium 1 869.45 0.9228 0.580
## + Barium 1 869.46 0.9115 0.595
## + Copper 1 869.46 0.9110 0.605
## + Nickel 1 869.55 0.8407 0.780
## + Moisture_Content 1 869.71 0.7150 0.880
## + Calcium 1 869.73 0.6989 0.880
## + Sulfur 1 869.67 0.7431 0.890
## + Organic_Matter 1 869.68 0.7370 0.905
## + Plot 0 868.62
## + total_precip_irr 0 868.62
m1_nifH_cca_chem <- cca(abund_tablenifH ~ Organic_Matter + Moisture_Content + Nitrate_Nitrite + Ammonia + Av_Phosphorus + Av_Potassium + pH + Barium + Calcium + Cobalt + Copper + Iron + Magnesium + Manganese + Nickel + Phosphorus + Sulfur + Zinc , meta3)
m0_nifH_cca_chem <- cca(abund_tablenifH ~ 1, meta3)
m1_nifH_cca_chem
## Call: cca(formula = abund_tablenifH ~ Organic_Matter +
## Moisture_Content + Nitrate_Nitrite + Ammonia + Av_Phosphorus +
## Av_Potassium + pH + Barium + Calcium + Cobalt + Copper + Iron +
## Magnesium + Manganese + Nickel + Phosphorus + Sulfur + Zinc, data
## = meta3)
##
## Inertia Proportion Rank
## Total 10.7659 1.0000
## Constrained 5.2066 0.4836 18
## Unconstrained 5.5593 0.5164 35
## Inertia is scaled Chi-square
##
## Eigenvalues for constrained axes:
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6 CCA7 CCA8 CCA9 CCA10
## 0.7916 0.6368 0.5430 0.5237 0.4736 0.4137 0.3403 0.2283 0.2008 0.1907
## CCA11 CCA12 CCA13 CCA14 CCA15 CCA16 CCA17 CCA18
## 0.1609 0.1430 0.1267 0.1085 0.0994 0.0838 0.0794 0.0625
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.4573 0.4358 0.3745 0.3252 0.2860 0.2803 0.2665 0.2449
## (Showed only 8 of all 35 unconstrained eigenvalues)
m0_nifH_cca_chem
## Call: cca(formula = abund_tablenifH ~ 1, data = meta3)
##
## Inertia Rank
## Total 10.77
## Unconstrained 10.77 53
## Inertia is scaled Chi-square
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.8189 0.6921 0.6401 0.5860 0.5584 0.5037 0.4818 0.4700
## (Showed only 8 of all 53 unconstrained eigenvalues)
model_nifH_cca_chem <-ordistep(m0_nifH_cca_chem, scope=formula(m1_nifH_cca_chem))
##
## Start: abund_tablenifH ~ 1
##
## Df AIC F Pr(>F)
## + Organic_Matter 1 874.75 3.5957 0.005 **
## + Nitrate_Nitrite 1 875.40 2.9283 0.005 **
## + Zinc 1 875.68 2.6463 0.005 **
## + Av_Potassium 1 875.71 2.6178 0.005 **
## + Nickel 1 875.88 2.4512 0.005 **
## + Cobalt 1 875.92 2.4074 0.005 **
## + Barium 1 875.93 2.3929 0.005 **
## + Iron 1 875.93 2.3923 0.005 **
## + Copper 1 875.99 2.3375 0.005 **
## + Magnesium 1 876.00 2.3264 0.005 **
## + Av_Phosphorus 1 876.01 2.3202 0.005 **
## + Phosphorus 1 876.07 2.2549 0.005 **
## + Moisture_Content 1 876.18 2.1464 0.005 **
## + Manganese 1 876.18 2.1462 0.005 **
## + Ammonia 1 876.47 1.8528 0.005 **
## + Calcium 1 876.57 1.7584 0.005 **
## + Sulfur 1 876.47 1.8514 0.010 **
## + pH 1 876.84 1.4897 0.030 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Organic_Matter
##
## Df AIC F Pr(>F)
## - Organic_Matter 1 876.36 3.5957 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Nitrate_Nitrite 1 873.56 3.1042 0.005 **
## + Zinc 1 873.97 2.6966 0.005 **
## + Barium 1 874.16 2.5125 0.005 **
## + Magnesium 1 874.20 2.4698 0.005 **
## + Cobalt 1 874.21 2.4614 0.005 **
## + Copper 1 874.22 2.4475 0.005 **
## + Av_Potassium 1 874.23 2.4431 0.005 **
## + Iron 1 874.24 2.4244 0.005 **
## + Av_Phosphorus 1 874.27 2.3965 0.005 **
## + Phosphorus 1 874.34 2.3325 0.005 **
## + Nickel 1 874.36 2.3098 0.005 **
## + Manganese 1 874.39 2.2789 0.005 **
## + Moisture_Content 1 874.48 2.1872 0.005 **
## + Sulfur 1 874.63 2.0475 0.005 **
## + Calcium 1 874.79 1.8862 0.005 **
## + Ammonia 1 874.88 1.7998 0.005 **
## + pH 1 875.14 1.5430 0.010 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite
##
## Df AIC F Pr(>F)
## - Nitrate_Nitrite 1 874.75 3.1042 0.005 **
## - Organic_Matter 1 875.40 3.7617 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Av_Potassium 1 872.79 2.6309 0.005 **
## + Magnesium 1 872.92 2.5090 0.005 **
## + Av_Phosphorus 1 872.98 2.4507 0.005 **
## + Phosphorus 1 872.99 2.4422 0.005 **
## + Moisture_Content 1 873.15 2.2870 0.005 **
## + Manganese 1 873.22 2.2192 0.005 **
## + Zinc 1 873.36 2.0786 0.005 **
## + Sulfur 1 873.39 2.0553 0.005 **
## + Copper 1 873.50 1.9447 0.005 **
## + Cobalt 1 873.54 1.9069 0.005 **
## + Calcium 1 873.55 1.8973 0.005 **
## + Barium 1 873.61 1.8369 0.005 **
## + Ammonia 1 873.65 1.8007 0.005 **
## + Iron 1 873.67 1.7782 0.005 **
## + Nickel 1 873.84 1.6164 0.005 **
## + pH 1 873.98 1.4887 0.020 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium
##
## Df AIC F Pr(>F)
## - Av_Potassium 1 873.56 2.6309 0.005 **
## - Nitrate_Nitrite 1 874.23 3.2820 0.005 **
## - Organic_Matter 1 874.51 3.5682 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Av_Phosphorus 1 872.25 2.3577 0.005 **
## + Magnesium 1 872.36 2.2606 0.005 **
## + Manganese 1 872.56 2.0639 0.005 **
## + Sulfur 1 872.61 2.0249 0.005 **
## + Calcium 1 872.63 2.0019 0.005 **
## + Ammonia 1 872.75 1.8850 0.005 **
## + Barium 1 872.80 1.8420 0.005 **
## + Nickel 1 872.90 1.7454 0.005 **
## + pH 1 872.93 1.7153 0.005 **
## + Copper 1 872.96 1.6867 0.005 **
## + Zinc 1 872.98 1.6722 0.005 **
## + Phosphorus 1 873.06 1.5945 0.005 **
## + Cobalt 1 872.89 1.7581 0.010 **
## + Moisture_Content 1 873.18 1.4875 0.015 *
## + Iron 1 873.16 1.4996 0.035 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium + Av_Phosphorus
##
## Df AIC F Pr(>F)
## - Av_Phosphorus 1 872.79 2.3577 0.005 **
## - Av_Potassium 1 872.98 2.5341 0.005 **
## - Nitrate_Nitrite 1 873.69 3.2208 0.005 **
## - Organic_Matter 1 874.15 3.6662 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Manganese 1 871.90 2.1402 0.005 **
## + Magnesium 1 872.03 2.0192 0.005 **
## + Sulfur 1 872.11 1.9489 0.005 **
## + Calcium 1 872.19 1.8727 0.005 **
## + Barium 1 872.29 1.7817 0.005 **
## + Ammonia 1 872.36 1.7145 0.005 **
## + Zinc 1 872.38 1.6967 0.005 **
## + Copper 1 872.44 1.6399 0.005 **
## + Phosphorus 1 872.48 1.6044 0.005 **
## + Cobalt 1 872.28 1.7883 0.010 **
## + Nickel 1 872.47 1.6161 0.010 **
## + Iron 1 872.54 1.5510 0.010 **
## + pH 1 872.58 1.5087 0.025 *
## + Moisture_Content 1 872.82 1.2883 0.080 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium + Av_Phosphorus + Manganese
##
## Df AIC F Pr(>F)
## - Manganese 1 872.25 2.1402 0.005 **
## - Av_Phosphorus 1 872.56 2.4287 0.005 **
## - Av_Potassium 1 872.74 2.5890 0.005 **
## - Nitrate_Nitrite 1 873.12 2.9469 0.005 **
## - Organic_Matter 1 873.85 3.6433 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Barium 1 871.44 2.1898 0.005 **
## + Sulfur 1 871.71 1.9435 0.005 **
## + Magnesium 1 871.76 1.9029 0.005 **
## + Zinc 1 871.79 1.8710 0.005 **
## + Cobalt 1 871.82 1.8413 0.005 **
## + Iron 1 871.91 1.7662 0.005 **
## + Copper 1 871.91 1.7627 0.005 **
## + Calcium 1 871.92 1.7571 0.005 **
## + Nickel 1 872.03 1.6558 0.005 **
## + Phosphorus 1 872.05 1.6338 0.005 **
## + pH 1 872.11 1.5817 0.005 **
## + Ammonia 1 872.20 1.5044 0.010 **
## + Moisture_Content 1 872.38 1.3407 0.050 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium + Av_Phosphorus + Manganese + Barium
##
## Df AIC F Pr(>F)
## - Barium 1 871.90 2.1898 0.005 **
## - Av_Phosphorus 1 872.08 2.3589 0.005 **
## - Nitrate_Nitrite 1 872.21 2.4735 0.005 **
## - Av_Potassium 1 872.29 2.5438 0.005 **
## - Manganese 1 872.29 2.5440 0.005 **
## - Organic_Matter 1 873.10 3.2928 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Cobalt 1 871.21 1.9353 0.005 **
## + Copper 1 871.30 1.8613 0.005 **
## + Sulfur 1 871.30 1.8589 0.005 **
## + Iron 1 871.31 1.8503 0.005 **
## + Magnesium 1 871.35 1.8125 0.005 **
## + Zinc 1 871.42 1.7495 0.005 **
## + Calcium 1 871.43 1.7405 0.005 **
## + Phosphorus 1 871.56 1.6326 0.005 **
## + Nickel 1 871.65 1.5524 0.005 **
## + pH 1 871.56 1.6330 0.010 **
## + Moisture_Content 1 871.84 1.3860 0.015 *
## + Ammonia 1 871.84 1.3831 0.025 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium + Av_Phosphorus + Manganese + Barium + Cobalt
##
## Df AIC F Pr(>F)
## - Cobalt 1 871.44 1.9353 0.005 **
## - Barium 1 871.82 2.2773 0.005 **
## - Av_Phosphorus 1 871.87 2.3225 0.005 **
## - Nitrate_Nitrite 1 871.91 2.3545 0.005 **
## - Manganese 1 871.92 2.3678 0.005 **
## - Av_Potassium 1 872.23 2.6409 0.005 **
## - Organic_Matter 1 873.08 3.4110 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Sulfur 1 870.86 2.0053 0.005 **
## + Calcium 1 870.96 1.9217 0.005 **
## + Nickel 1 871.19 1.7216 0.005 **
## + pH 1 871.48 1.4699 0.010 **
## + Phosphorus 1 871.57 1.3904 0.025 *
## + Moisture_Content 1 871.58 1.3814 0.025 *
## + Ammonia 1 871.61 1.3613 0.040 *
## + Zinc 1 871.68 1.2995 0.060 .
## + Magnesium 1 871.71 1.2731 0.070 .
## + Copper 1 871.87 1.1341 0.215
## + Iron 1 872.01 1.0150 0.530
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium + Av_Phosphorus + Manganese + Barium + Cobalt + Sulfur
##
## Df AIC F Pr(>F)
## - Sulfur 1 871.21 2.0053 0.005 **
## - Cobalt 1 871.30 2.0804 0.005 **
## - Barium 1 871.39 2.1625 0.005 **
## - Av_Phosphorus 1 871.44 2.1990 0.005 **
## - Manganese 1 871.57 2.3147 0.005 **
## - Nitrate_Nitrite 1 871.67 2.4073 0.005 **
## - Av_Potassium 1 872.03 2.7182 0.005 **
## - Organic_Matter 1 872.99 3.5736 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + Magnesium 1 871.12 1.4435 0.010 **
## + Nickel 1 871.12 1.4399 0.040 *
## + pH 1 871.36 1.2417 0.090 .
## + Copper 1 871.38 1.2255 0.125
## + Phosphorus 1 871.35 1.2508 0.130
## + Zinc 1 871.41 1.2013 0.145
## + Ammonia 1 871.53 1.0985 0.305
## + Moisture_Content 1 871.63 1.0163 0.485
## + Iron 1 871.66 0.9862 0.520
## + Calcium 1 871.93 0.7671 0.835
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium + Av_Phosphorus + Manganese + Barium + Cobalt + Sulfur + Magnesium
##
## Df AIC F Pr(>F)
## - Cobalt 1 870.55 1.1831 0.225
## - Magnesium 1 870.86 1.4435 0.030 *
## - Av_Phosphorus 1 870.94 1.5081 0.005 **
## - Barium 1 871.20 1.7285 0.005 **
## - Sulfur 1 871.71 2.1626 0.005 **
## - Manganese 1 871.72 2.1755 0.005 **
## - Nitrate_Nitrite 1 871.80 2.2421 0.005 **
## - Organic_Matter 1 871.94 2.3606 0.005 **
## - Av_Potassium 1 872.38 2.7380 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium + Av_Phosphorus + Manganese + Barium + Sulfur + Magnesium
##
## Df AIC F Pr(>F)
## + Phosphorus 1 870.39 1.7956 0.005 **
## + Iron 1 870.99 1.2859 0.155
## + Ammonia 1 871.17 1.1384 0.210
## + pH 1 871.19 1.1215 0.230
## + Cobalt 1 871.12 1.1831 0.245
## + Nickel 1 871.20 1.1160 0.310
## + Zinc 1 871.44 0.9148 0.650
## + Moisture_Content 1 871.48 0.8821 0.675
## + Copper 1 871.45 0.9088 0.710
## + Calcium 1 871.50 0.8643 0.740
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium + Av_Phosphorus + Manganese + Barium + Sulfur + Magnesium + Phosphorus
##
## Df AIC F Pr(>F)
## - Av_Phosphorus 1 870.15 1.4610 0.025 *
## - Av_Potassium 1 870.27 1.5617 0.010 **
## - Phosphorus 1 870.55 1.7956 0.005 **
## - Magnesium 1 870.72 1.9382 0.005 **
## - Nitrate_Nitrite 1 870.93 2.1213 0.005 **
## - Sulfur 1 871.03 2.2060 0.005 **
## - Manganese 1 871.07 2.2387 0.005 **
## - Barium 1 871.10 2.2666 0.005 **
## - Organic_Matter 1 871.48 2.5897 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Df AIC F Pr(>F)
## + pH 1 870.92 1.1871 0.180
## + Ammonia 1 870.94 1.1697 0.215
## + Nickel 1 871.03 1.0927 0.275
## + Cobalt 1 871.04 1.0916 0.275
## + Iron 1 871.03 1.0935 0.315
## + Zinc 1 871.13 1.0140 0.445
## + Calcium 1 871.23 0.9362 0.510
## + Moisture_Content 1 871.28 0.8899 0.735
## + Copper 1 871.33 0.8507 0.775
model_nifH_cca_chem$anova
## Df AIC F Pr(>F)
## + Organic_Matter 1 874.75 3.5957 0.005 **
## + Nitrate_Nitrite 1 873.56 3.1042 0.005 **
## + Av_Potassium 1 872.79 2.6309 0.005 **
## + Av_Phosphorus 1 872.25 2.3577 0.005 **
## + Manganese 1 871.90 2.1402 0.005 **
## + Barium 1 871.44 2.1898 0.005 **
## + Cobalt 1 871.21 1.9353 0.005 **
## + Sulfur 1 870.86 2.0053 0.005 **
## + Magnesium 1 871.12 1.4435 0.010 **
## - Cobalt 1 870.55 1.1831 0.225
## + Phosphorus 1 870.39 1.7956 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1_nifH_rda <- rda(abund_tablenifH ~ ., meta3)
m0_nifH_rda <- rda(abund_tablenifH ~ 1, meta3)
m1_nifH_rda
## Call: rda(formula = abund_tablenifH ~ Site + Pea_variety + Plot +
## season_precip + irrgation + total_precip_irr + Tillage + prev_crop
## + grain_yield + Organic_Matter + Moisture_Content +
## Nitrate_Nitrite + Ammonia + Av_Phosphorus + Av_Potassium + pH +
## Barium + Calcium + Cobalt + Copper + Iron + Magnesium + Manganese
## + Nickel + Phosphorus + Sulfur + Zinc, data = meta3)
##
## Inertia Proportion Rank
## Total 1.258e+11 1.000e+00
## Constrained 9.747e+10 7.745e-01 33
## Unconstrained 2.837e+10 2.255e-01 20
## Inertia is variance
## Some constraints were aliased because they were collinear (redundant)
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6
## 16353665417 15356453441 12209872418 10446884520 7322044854 5776558730
## RDA7 RDA8 RDA9 RDA10 RDA11 RDA12
## 4833284412 4476474731 3654563170 3048908811 2374942944 1890413915
## RDA13 RDA14 RDA15 RDA16 RDA17 RDA18
## 1488069677 1072841995 989712375 941967373 836809708 626843871
## RDA19 RDA20 RDA21 RDA22 RDA23 RDA24
## 573856607 508536865 403729752 376714333 325054973 282465322
## RDA25 RDA26 RDA27 RDA28 RDA29 RDA30
## 238159099 207287067 185699351 163253845 127967489 116333192
## RDA31 RDA32 RDA33
## 95605085 84351951 75938358
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6
## 5599447125 4342938458 3880332013 3368259818 2027561149 1483067201
## PC7 PC8
## 1310478544 1156592863
## (Showed only 8 of all 20 unconstrained eigenvalues)
m0_nifH_rda
## Call: rda(formula = abund_tablenifH ~ 1, data = meta3)
##
## Inertia Rank
## Total 1.258e+11
## Unconstrained 1.258e+11 53
## Inertia is variance
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6
## 18388165460 17017988044 13929669751 13228540294 10298920752 7474983345
## PC7 PC8
## 6353906098 5386444205
## (Showed only 8 of all 53 unconstrained eigenvalues)
plot(m1_nifH_rda)
model_rda_nifH <-ordiR2step(m0_nifH_rda, scope=formula(m1_nifH_rda))
## Step: R2.adj= 0
## Call: abund_tablenifH ~ 1
##
## R2.adjusted
## <All variables> 0.40253072
## + Site 0.39284983
## + prev_crop 0.14946826
## + Tillage 0.09118753
## + total_precip_irr 0.06245835
## + Av_Phosphorus 0.06052058
## + Phosphorus 0.05195778
## + Av_Potassium 0.05112278
## + Magnesium 0.04894792
## + irrgation 0.04678398
## + Plot 0.04077308
## + Moisture_Content 0.04047656
## + season_precip 0.03754349
## + Barium 0.03466802
## + Manganese 0.03434338
## + Zinc 0.03212237
## + Organic_Matter 0.03177594
## + Ammonia 0.02812448
## + Nitrate_Nitrite 0.02759831
## + Nickel 0.02460517
## + Copper 0.02201680
## + Sulfur 0.02175055
## + Calcium 0.01980531
## + grain_yield 0.01948985
## + Iron 0.01591920
## + pH 0.01585647
## + Cobalt 0.01268035
## <none> 0.00000000
## + Pea_variety -0.02472576
##
## Df AIC F Pr(>F)
## + Site 8 1361.4 5.2866 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.3928498
## Call: abund_tablenifH ~ Site
##
## R2.adjusted
## + Pea_variety 0.4053417
## <All variables> 0.4025307
## + Ammonia 0.3971916
## + Nitrate_Nitrite 0.3945793
## <none> 0.3928498
## + Plot 0.3928498
## + season_precip 0.3928498
## + irrgation 0.3928498
## + total_precip_irr 0.3928498
## + Tillage 0.3928498
## + prev_crop 0.3928498
## + Av_Potassium 0.3920019
## + Organic_Matter 0.3905573
## + Av_Phosphorus 0.3888507
## + Sulfur 0.3888353
## + pH 0.3883561
## + Calcium 0.3881788
## + Moisture_Content 0.3878418
## + Manganese 0.3866273
## + grain_yield 0.3860617
## + Magnesium 0.3853474
## + Nickel 0.3852829
## + Cobalt 0.3846202
## + Copper 0.3845997
## + Iron 0.3842913
## + Barium 0.3842156
## + Phosphorus 0.3841708
## + Zinc 0.3838492
m1_nifH_rda_chem <- rda(abund_tablenifH ~ Organic_Matter + Moisture_Content + Nitrate_Nitrite + Ammonia + Av_Phosphorus + Av_Potassium + pH + Barium + Calcium + Cobalt + Copper + Iron + Magnesium + Manganese + Nickel + Phosphorus + Sulfur + Zinc , meta3)
m0_nifH_rda_chem <- rda(abund_tablenifH ~ 1, meta3)
m1_nifH_rda_chem
## Call: rda(formula = abund_tablenifH ~ Organic_Matter +
## Moisture_Content + Nitrate_Nitrite + Ammonia + Av_Phosphorus +
## Av_Potassium + pH + Barium + Calcium + Cobalt + Copper + Iron +
## Magnesium + Manganese + Nickel + Phosphorus + Sulfur + Zinc, data
## = meta3)
##
## Inertia Proportion Rank
## Total 1.258e+11 1.000e+00
## Constrained 6.893e+10 5.478e-01 18
## Unconstrained 5.690e+10 4.522e-01 35
## Inertia is variance
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6
## 14804112412 13108358521 10437172065 8136887370 5686929136 3911450523
## RDA7 RDA8 RDA9 RDA10 RDA11 RDA12
## 3487763724 2618476736 1896564759 1302627898 875978909 653782995
## RDA13 RDA14 RDA15 RDA16 RDA17 RDA18
## 566194235 425743412 312506701 290536181 223053545 193225475
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6
## 11018590525 7155476437 5939185648 5422891779 4017288576 3016129382
## PC7 PC8
## 2886770730 2071617511
## (Showed only 8 of all 35 unconstrained eigenvalues)
m0_nifH_rda_chem
## Call: rda(formula = abund_tablenifH ~ 1, data = meta3)
##
## Inertia Rank
## Total 1.258e+11
## Unconstrained 1.258e+11 53
## Inertia is variance
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6
## 18388165460 17017988044 13929669751 13228540294 10298920752 7474983345
## PC7 PC8
## 6353906098 5386444205
## (Showed only 8 of all 53 unconstrained eigenvalues)
plot(m1_nifH_rda_chem)
model_rda_chem_nifH <-ordiR2step(m0_nifH_rda_chem, scope=formula(m1_nifH_rda_chem))
## Step: R2.adj= 0
## Call: abund_tablenifH ~ 1
##
## R2.adjusted
## <All variables> 0.31521855
## + Av_Phosphorus 0.06052058
## + Phosphorus 0.05195778
## + Av_Potassium 0.05112278
## + Magnesium 0.04894792
## + Moisture_Content 0.04047656
## + Barium 0.03466802
## + Manganese 0.03434338
## + Zinc 0.03212237
## + Organic_Matter 0.03177594
## + Ammonia 0.02812448
## + Nitrate_Nitrite 0.02759831
## + Nickel 0.02460517
## + Copper 0.02201680
## + Sulfur 0.02175055
## + Calcium 0.01980531
## + Iron 0.01591920
## + pH 0.01585647
## + Cobalt 0.01268035
## <none> 0.00000000
##
## Df AIC F Pr(>F)
## + Av_Phosphorus 1 1378.7 4.4142 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.06052058
## Call: abund_tablenifH ~ Av_Phosphorus
##
## R2.adjusted
## <All variables> 0.31521855
## + Av_Potassium 0.10253413
## + Phosphorus 0.09795318
## + Zinc 0.09505177
## + Manganese 0.09455271
## + Organic_Matter 0.09324882
## + Ammonia 0.09270795
## + Nickel 0.08781964
## + Magnesium 0.08490752
## + Barium 0.08424954
## + Moisture_Content 0.08392402
## + Nitrate_Nitrite 0.08340692
## + Sulfur 0.08244671
## + Copper 0.08171114
## + Iron 0.08134433
## + Calcium 0.07884137
## + Cobalt 0.07764768
## + pH 0.07328876
## <none> 0.06052058
##
## Df AIC F Pr(>F)
## + Av_Potassium 1 1377.2 3.4343 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.1025341
## Call: abund_tablenifH ~ Av_Phosphorus + Av_Potassium
##
## R2.adjusted
## <All variables> 0.3152185
## + Organic_Matter 0.1426010
## + Zinc 0.1346765
## + Manganese 0.1346442
## + Sulfur 0.1321684
## + Nickel 0.1314652
## + Phosphorus 0.1297865
## + Calcium 0.1276961
## + Magnesium 0.1269648
## + Nitrate_Nitrite 0.1268511
## + Ammonia 0.1264398
## + Barium 0.1263189
## + Copper 0.1232790
## + Moisture_Content 0.1227996
## + pH 0.1220798
## + Iron 0.1192709
## + Cobalt 0.1140726
## <none> 0.1025341
##
## Df AIC F Pr(>F)
## + Organic_Matter 1 1375.7 3.3833 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.142601
## Call: abund_tablenifH ~ Av_Phosphorus + Av_Potassium + Organic_Matter
##
## R2.adjusted
## <All variables> 0.3152185
## + Barium 0.1806009
## + Manganese 0.1778531
## + Zinc 0.1739765
## + Ammonia 0.1728053
## + Sulfur 0.1698627
## + Calcium 0.1691275
## + Magnesium 0.1685205
## + Nitrate_Nitrite 0.1637934
## + Nickel 0.1637017
## + Copper 0.1636724
## + pH 0.1635979
## + Phosphorus 0.1610160
## + Iron 0.1572052
## + Moisture_Content 0.1536428
## + Cobalt 0.1526631
## <none> 0.1426010
##
## Df AIC F Pr(>F)
## + Barium 1 1374.1 3.3188 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.1806009
## Call: abund_tablenifH ~ Av_Phosphorus + Av_Potassium + Organic_Matter + Barium
##
## R2.adjusted
## <All variables> 0.3152185
## + Manganese 0.2425367
## + Sulfur 0.2152786
## + Calcium 0.2146288
## + Magnesium 0.2075492
## + Cobalt 0.2059481
## + Iron 0.2045766
## + Zinc 0.2041698
## + Nickel 0.2040281
## + pH 0.2023023
## + Ammonia 0.2016233
## + Copper 0.2015660
## + Phosphorus 0.2011897
## + Moisture_Content 0.1920917
## + Nitrate_Nitrite 0.1880401
## <none> 0.1806009
##
## Df AIC F Pr(>F)
## + Manganese 1 1370.8 5.0066 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.2425367
## Call: abund_tablenifH ~ Av_Phosphorus + Av_Potassium + Organic_Matter + Barium + Manganese
##
## R2.adjusted
## <All variables> 0.3152185
## + Iron 0.2741671
## + Magnesium 0.2715470
## + Zinc 0.2703388
## + Sulfur 0.2688961
## + Copper 0.2683106
## + Cobalt 0.2675120
## + pH 0.2661270
## + Calcium 0.2650136
## + Phosphorus 0.2645607
## + Nickel 0.2626842
## + Moisture_Content 0.2555896
## + Ammonia 0.2539120
## + Nitrate_Nitrite 0.2511329
## <none> 0.2425367
##
## Df AIC F Pr(>F)
## + Iron 1 1369.3 3.0917 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.2741671
## Call: abund_tablenifH ~ Av_Phosphorus + Av_Potassium + Organic_Matter + Barium + Manganese + Iron
##
## R2.adjusted
## <All variables> 0.3152185
## + Sulfur 0.3077588
## + Calcium 0.3034844
## + pH 0.2895106
## + Phosphorus 0.2873660
## + Moisture_Content 0.2869075
## + Ammonia 0.2867278
## + Nickel 0.2850464
## + Zinc 0.2833410
## + Magnesium 0.2802296
## + Nitrate_Nitrite 0.2787387
## + Cobalt 0.2743731
## <none> 0.2741671
## + Copper 0.2741048
##
## Df AIC F Pr(>F)
## + Sulfur 1 1367.6 3.2807 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.3077588
## Call: abund_tablenifH ~ Av_Phosphorus + Av_Potassium + Organic_Matter + Barium + Manganese + Iron + Sulfur
##
## R2.adjusted
## + Magnesium 0.3193868
## + Zinc 0.3191464
## <All variables> 0.3152185
## + Ammonia 0.3141646
## + pH 0.3132794
## + Phosphorus 0.3131892
## + Nitrate_Nitrite 0.3129262
## + Copper 0.3093376
## + Moisture_Content 0.3090885
## + Nickel 0.3081999
## <none> 0.3077588
## + Cobalt 0.3049041
## + Calcium 0.3007448
model_rda_chem_nifH$anova
## R2.adj Df AIC F Pr(>F)
## + Av_Phosphorus 0.060521 1 1378.7 4.4142 0.002 **
## + Av_Potassium 0.102534 1 1377.2 3.4343 0.002 **
## + Organic_Matter 0.142601 1 1375.7 3.3833 0.002 **
## + Barium 0.180601 1 1374.1 3.3188 0.002 **
## + Manganese 0.242537 1 1370.8 5.0066 0.002 **
## + Iron 0.274167 1 1369.3 3.0917 0.002 **
## + Sulfur 0.307759 1 1367.6 3.2807 0.002 **
## <All variables> 0.315219
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1_nifH_cap_chem<- capscale(abund_tablenifH ~ Organic_Matter + Moisture_Content + Nitrate_Nitrite + Ammonia + Av_Phosphorus + Av_Potassium + pH + Barium + Calcium + Cobalt + Copper + Iron + Magnesium + Manganese + Nickel + Phosphorus + Sulfur + Zinc , data = meta3, distance = "bray")
m0_nifH_cap_chem<- capscale(abund_distnifH ~ 1, meta3)
plot(m1_nifH_cap_chem)
model_cap_chem_nifH <-ordiR2step(m0_nifH_cap_chem, scope=formula(m1_nifH_cap_chem))
## Step: R2.adj= 0
## Call: abund_distnifH ~ 1
##
## R2.adjusted
## <All variables> 0.42412042
## + Organic_Matter 0.08664823
## + Av_Potassium 0.06424058
## + Av_Phosphorus 0.05197708
## + Zinc 0.04843832
## + Nitrate_Nitrite 0.04757351
## + Nickel 0.04703479
## + Barium 0.04533845
## + Iron 0.04135991
## + Magnesium 0.04111715
## + Cobalt 0.03937854
## + Phosphorus 0.03739460
## + Manganese 0.03679948
## + Copper 0.03582041
## + Ammonia 0.03314279
## + Moisture_Content 0.02637968
## + Sulfur 0.02145207
## + pH 0.02106463
## + Calcium 0.02096101
## <none> 0.00000000
##
## Df AIC F Pr(>F)
## + Organic_Matter 1 159.43 5.8777 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.08664823
## Call: abund_distnifH ~ Organic_Matter
##
## R2.adjusted
## <All variables> 0.42412042
## + Barium 0.14270722
## + Av_Potassium 0.14238286
## + Av_Phosphorus 0.13935754
## + Nitrate_Nitrite 0.13484559
## + Zinc 0.13050213
## + Magnesium 0.12895812
## + Phosphorus 0.12745386
## + Manganese 0.12724465
## + Copper 0.12206419
## + Iron 0.12132345
## + Cobalt 0.12126082
## + Ammonia 0.11844446
## + Nickel 0.11807834
## + Sulfur 0.11669887
## + Moisture_Content 0.11655745
## + Calcium 0.11156606
## + pH 0.10979577
## <none> 0.08664823
##
## Df AIC F Pr(>F)
## + Barium 1 157.08 4.2814 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.1427072
## Call: abund_distnifH ~ Organic_Matter + Barium
##
## R2.adjusted
## <All variables> 0.4241204
## + Av_Potassium 0.2008529
## + Av_Phosphorus 0.1900137
## + Manganese 0.1880846
## + Phosphorus 0.1853480
## + Magnesium 0.1842902
## + Cobalt 0.1825912
## + Copper 0.1820141
## + Zinc 0.1809508
## + Iron 0.1802010
## + Ammonia 0.1754906
## + Moisture_Content 0.1750410
## + Nickel 0.1717581
## + Sulfur 0.1707810
## + Calcium 0.1678310
## + pH 0.1675754
## + Nitrate_Nitrite 0.1669577
## <none> 0.1427072
##
## Df AIC F Pr(>F)
## + Av_Potassium 1 154.35 4.5718 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.2008529
## Call: abund_distnifH ~ Organic_Matter + Barium + Av_Potassium
##
## R2.adjusted
## <All variables> 0.4241204
## + Manganese 0.2571847
## + Av_Phosphorus 0.2465866
## + Cobalt 0.2426324
## + Copper 0.2399983
## + Magnesium 0.2390988
## + Iron 0.2386871
## + Ammonia 0.2363363
## + Zinc 0.2360117
## + Nickel 0.2320877
## + pH 0.2311367
## + Sulfur 0.2301444
## + Calcium 0.2287892
## + Phosphorus 0.2236127
## + Nitrate_Nitrite 0.2227016
## + Moisture_Content 0.2141867
## <none> 0.2008529
##
## Df AIC F Pr(>F)
## + Manganese 1 151.47 4.637 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.2571847
## Call: abund_distnifH ~ Organic_Matter + Barium + Av_Potassium + Manganese
##
## R2.adjusted
## <All variables> 0.4241204
## + Av_Phosphorus 0.3048552
## + Copper 0.2999316
## + Magnesium 0.2990454
## + Iron 0.2975378
## + Zinc 0.2954468
## + Cobalt 0.2944697
## + pH 0.2893062
## + Nickel 0.2881294
## + Ammonia 0.2861456
## + Sulfur 0.2846872
## + Calcium 0.2815357
## + Phosphorus 0.2812081
## + Nitrate_Nitrite 0.2764049
## + Moisture_Content 0.2720061
## <none> 0.2571847
##
## Df AIC F Pr(>F)
## + Av_Phosphorus 1 148.93 4.2071 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.3048552
## Call: abund_distnifH ~ Organic_Matter + Barium + Av_Potassium + Manganese + Av_Phosphorus
##
## R2.adjusted
## <All variables> 0.4241204
## + Iron 0.3478017
## + Cobalt 0.3450564
## + Copper 0.3442042
## + Magnesium 0.3418955
## + Zinc 0.3394628
## + pH 0.3331523
## + Sulfur 0.3315043
## + Phosphorus 0.3300620
## + Nickel 0.3283710
## + Calcium 0.3275764
## + Nitrate_Nitrite 0.3265021
## + Moisture_Content 0.3154696
## + Ammonia 0.3126530
## <none> 0.3048552
##
## Df AIC F Pr(>F)
## + Iron 1 146.52 4.0021 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.3478017
## Call: abund_distnifH ~ Organic_Matter + Barium + Av_Potassium + Manganese + Av_Phosphorus + Iron
##
## R2.adjusted
## <All variables> 0.4241204
## + Sulfur 0.3885315
## + Calcium 0.3839583
## + Nickel 0.3693593
## + pH 0.3689149
## + Nitrate_Nitrite 0.3651795
## + Zinc 0.3638021
## + Phosphorus 0.3632509
## + Moisture_Content 0.3622189
## + Magnesium 0.3594648
## + Ammonia 0.3562554
## + Copper 0.3550172
## + Cobalt 0.3480682
## <none> 0.3478017
##
## Df AIC F Pr(>F)
## + Sulfur 1 144.06 3.9597 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.3885315
## Call: abund_distnifH ~ Organic_Matter + Barium + Av_Potassium + Manganese + Av_Phosphorus + Iron + Sulfur
##
## R2.adjusted
## <All variables> 0.4241204
## + Magnesium 0.4112957
## + Zinc 0.4085774
## + Nitrate_Nitrite 0.4069902
## + pH 0.4051098
## + Nickel 0.4019619
## + Phosphorus 0.4015875
## + Copper 0.3991158
## + Moisture_Content 0.3920267
## + Ammonia 0.3894283
## <none> 0.3885315
## + Calcium 0.3869749
## + Cobalt 0.3831228
##
## Df AIC F Pr(>F)
## + Magnesium 1 142.96 2.6569 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.4112957
## Call: abund_distnifH ~ Organic_Matter + Barium + Av_Potassium + Manganese + Av_Phosphorus + Iron + Sulfur + Magnesium
##
## R2.adjusted
## + Phosphorus 0.4318158
## + Nitrate_Nitrite 0.4268994
## <All variables> 0.4241204
## + Nickel 0.4189326
## + pH 0.4174772
## + Zinc 0.4166306
## + Ammonia 0.4127613
## <none> 0.4112957
## + Moisture_Content 0.4108080
## + Calcium 0.4093634
## + Copper 0.4068650
## + Cobalt 0.4064199
“http://deneflab.github.io/MicrobeMiseq/demos/mothur_2_phyloseq.html#constrained_ordinations”
# CCA ordinate
cca_ord_nifH <- ordinate(
physeq = physeq_nifH_ord,
method = "CCA",
distance = abund_distnifH,
formula = ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium +
Av_Phosphorus + Manganese + Barium + Sulfur + Magnesium )
# CCA plot
cca_plot_nifH <- plot_ordination(
physeq = physeq_nifH_ord,
ordination = cca_ord_nifH,
color = "Site",
axes = c(1,2)) +
geom_point(aes(colour = Site), size = 3) +
scale_color_manual(values = farm_col_paired)
# Now add the environmental variables as arrows
arrowmat_nifH_cca <- vegan::scores(cca_ord_nifH, display = "bp")
#Get appropiate scalling multipler
mul<-vegan::ordiArrowMul(arrowmat_nifH_cca)
#Multiply biplot by scaling multiplier
arrowmat_nifH_cca_scale<-arrowmat_nifH_cca*3
# Add labels, make a data.frame
arrowdf_nifH_cca <- data.frame(labels = rownames(arrowmat_nifH_cca_scale), arrowmat_nifH_cca_scale)
# Define the arrow aesthetic mapping
arrow_map <- aes(xend = CCA1,
yend = CCA2,
x = 0,
y = 0,
shape = NULL,
color = NULL,
label = labels)
label_map <- aes(x = 1.3* CCA1,
y = 1.3 * CCA2,
shape = NULL,
color = NULL,
label = labels)
arrowhead = arrow(length = unit(0.02, "npc"))
# Make a new graphic
cca_plot_nifH +
geom_segment(
mapping = arrow_map,
size = .5,
data = arrowdf_nifH_cca,
color = "gray",
arrow = arrowhead
) +
geom_text(
mapping = label_map,
size = 4,
data = arrowdf_nifH_cca,
show.legend = FALSE
)+
ggtitle("CCA plot constrained ordination of nifH with selected Chemistry Model")+
theme_bw()
## Warning: Ignoring unknown aesthetics: label
# RDA ordinate
rda_ord_nifH <- ordinate(
physeq = physeq_nifH_ord,
method = "RDA",
distance = abund_distnifH,
formula = ~ Av_Phosphorus + Sulfur + Magnesium + Zinc + Av_Potassium + Manganese )
# RDA plot
rda_plot_nifH <- plot_ordination(
physeq = physeq_nifH_ord,
ordination = rda_ord_nifH,
color = "Site",
axes = c(1,2)) +
geom_point(aes(colour = Site), size = 3) +
scale_color_manual(values = farm_col_paired)
# Now add the environmental variables as arrows
arrowmat_nifH_rda <- vegan::scores(rda_ord_nifH, display = "bp")
#Get appropiate scalling multipler
mul<-vegan::ordiArrowMul(arrowmat_nifH_rda)
#Multiply biplot by scaling multiplier
arrowmat_nifH_rda_scale<-arrowmat_nifH_rda*700
# Add labels, make a data.frame
arrowdf_nifH_rda <- data.frame(labels = rownames(arrowmat_nifH_rda_scale), arrowmat_nifH_rda_scale)
# Define the arrow aesthetic mapping
arrow_map <- aes(xend = RDA1,
yend = RDA2,
x = 0,
y = 0,
shape = NULL,
color = NULL,
label = labels)
label_map <- aes(x = 1.3* RDA1,
y = 1.3 * RDA2,
shape = NULL,
color = NULL,
label = labels)
arrowhead = arrow(length = unit(0.02, "npc"))
# Make a new graphic
rda_plot_nifH +
geom_segment(
mapping = arrow_map,
size = .5,
data = arrowdf_nifH_rda,
color = "gray",
arrow = arrowhead
) +
geom_text(
mapping = label_map,
size = 4,
data = arrowdf_nifH_rda,
show.legend = FALSE
)+
ggtitle("RDA plot constrained ordination of nifH with selected Chemistry Model")+
theme_bw()
## Warning: Ignoring unknown aesthetics: label
# CAP ordinate
cap_ord_nifH <- ordinate(
physeq = physeq_nifH_ord,
method = "CAP",
distance = abund_distnifH,
formula = ~ Organic_Matter + Barium + Av_Potassium +
Manganese + Av_Phosphorus + Iron + Sulfur + Magnesium )
# CAP plot
cap_plot_nifH <- plot_ordination(
physeq = physeq_nifH_ord,
ordination = cap_ord_nifH,
color = "Site",
axes = c(1,2)) +
geom_point(aes(colour = Site), size = 3) +
scale_color_manual(values = farm_col_paired)
# Now add the environmental variables as arrows
arrowmat_nifH_cap <- vegan::scores(cap_ord_nifH, display = "bp")
#Get appropiate scalling multipler
mul<-vegan::ordiArrowMul(arrowmat_nifH_cap)
#Multiply biplot by scaling multiplier
arrowmat_nifH_cap_scale<-arrowmat_nifH_cap*1.9
# Add labels, make a data.frame
arrowdf_nifH_cap <- data.frame(labels = rownames(arrowmat_nifH_cap_scale), arrowmat_nifH_cap_scale)
# Define the arrow aesthetic mapping
arrow_map <- aes(xend = CAP1,
yend = CAP2,
x = 0,
y = 0,
shape = NULL,
color = NULL,
label = labels)
label_map <- aes(x = 1.3 * CAP1,
y = 1.3 * CAP2,
shape = NULL,
color = NULL,
label = labels)
arrowhead = arrow(length = unit(0.02, "npc"))
# Make a new graphic
cap_plot_nifH +
geom_segment(
mapping = arrow_map,
size = .5,
data = arrowdf_nifH_cap,
color = "gray",
arrow = arrowhead
) +
geom_text(
mapping = label_map,
size = 4,
data = arrowdf_nifH_cap,
show.legend = FALSE
)+
ggtitle("CAP Plot Constrained Ordination of nifH with Selected Model")+
theme_bw()
## Warning: Ignoring unknown aesthetics: label